MODULE DISTANCE
!this module solves the optimal consumption by directly solving the maximization problem
!by Simplex method.
  USE GLOBVAR
  USE INTERPOL
  USE NRUTIL, only : NR_BIG2, NR_SMALL, NR_SMALL2, ARRAY2D, NR_EPS
  use comf_trans, only : real2ab_exp, ab2real_exp
  use STATISTICS, only : Mean, Standard_Deviation, Correlation
  use probability, only : CDF_Normal
  use ROUWENHORST
  use COMFUNCS
  use INIT
  USE CMAXEE
  USE SIMPLEXC

  use MPI  !For Fortran 90, working with Intel version of openmpi.
  IMPLICIT NONE
  !include 'mpif.h'  !place this after IMPLICIT NONE, otherwise debug reports error. 
  
    CONTAINS

  !==========================================================================  
  !distance to be minimized, given the set of parameters. 
  !this is for subplx; input parameters already transformed to whole R
  !called in MASTER    
  !==========================================================================
  real(8) function SMM_distance_nn(nn,theta_)
    implicit none
    integer, intent(in) :: nn
    REAL(8), INTENT(IN) :: theta_(nn)
    integer :: i
    real(8) :: theta_act_(N_PAR)
    
    if (gcL_searchptonly) then
        theta_act_ = gr_Vpar_orig
        do i = 1, N_plpt
            theta_act_(N_PAR-N_plpt+i) = theta_(i)    !26,27,28
        enddo !do i=1
        SMM_distance_nn = SMM_distance(theta_act_)
    elseif (gcL_searchSSAonly) then
        theta_act_ = gr_Vpar_orig
        do i = 1, N_ssa_u
            theta_act_(22+i) = theta_(i)    !23,24,25
        enddo !do i=1
        SMM_distance_nn = SMM_distance(theta_act_)    
    elseif (gcI_robust == 21) then !no H depreciation at all
        theta_act_(1) = theta_(1)
        theta_act_(2) = 0.0d0  !redundant, as it will be set again in SMM_distance
        do i = 2, N_PAR-1
            theta_act_(i+1) = theta_(i)
        enddo
        SMM_distance_nn = SMM_distance(theta_act_)
    elseif (gcI_robust==12 .OR. gcI_robust==13) then !exo/lbd Hone
        theta_act_(1:2) = theta_(1:2)
        theta_act_(3) = 0.0d0  !alpha_I, set as 0.0d0
        do i = 3, N_PAR-1
            theta_act_(i+1) = theta_(i)
        enddo
        SMM_distance_nn = SMM_distance(theta_act_)
    else
        SMM_distance_nn = SMM_distance(theta_)
    endif
    return
  end function SMM_distance_nn
  
  subroutine sub_SMM_distance(theta_,val)
    implicit none
    REAL(8), INTENT(IN), dimension(:) :: theta_
    REAL(8), INTENT(OUT) :: val
    integer :: i
    real(8) :: theta_act_(N_PAR)
    
    if (gcL_searchptonly) then
        theta_act_ = gr_Vpar_orig
        do i = 1, N_plpt
            theta_act_(N_PAR-N_plpt+i) = theta_(i)    !26,27,28
        enddo !do i=1
        val = SMM_distance(theta_act_)
    elseif (gcL_searchSSAonly) then
        theta_act_ = gr_Vpar_orig
        do i = 1, N_ssa_u
            theta_act_(22+i) = theta_(i)    !23,24,25
        enddo !do i=1
        val = SMM_distance(theta_act_)    
    elseif (gcI_robust == 21) then !no H depreciation at all
        theta_act_(1) = theta_(1)
        theta_act_(2) = 0.0d0  !redundant, as it will be set again in SMM_distance
        do i = 2, N_PAR-1
            theta_act_(i+1) = theta_(i)
        enddo
        val = SMM_distance(theta_act_)
    elseif (gcI_robust==12 .OR. gcI_robust==13) then !exo/lbd Hone
        theta_act_(1:2) = theta_(1:2)
        theta_act_(3) = 0.0d0  !alpha_I, set as 0.0d0
        do i = 3, N_PAR-1
            theta_act_(i+1) = theta_(i)
        enddo
        val = SMM_distance(theta_act_)
    else
        val = SMM_distance(theta_)
    endif
    return
  end subroutine
  
  function SMM_distance_diff(theta_)
    implicit none
    REAL(8), INTENT(IN) :: theta_(:)
    real(8) :: SMM_distance_diff(size(theta_(:),1))
    real(8) :: x(size(theta_(:),1)), r0, r1
    integer :: i
    r0 = SMM_distance(theta_)
    do i = 1, N_PAR
       x = theta_
       x(i) = (1.0d0+1.0d-3) * theta_(i)
       r1 = SMM_distance(x)
       SMM_distance_diff(i) = (r1-r0)/(x(i)-theta_(i))
    enddo
    return
  end function SMM_distance_diff
  
  !==========================================================================  
  !called in MASTER    
  !==========================================================================
  real(8) function SMM_distance(theta_)
    implicit none
    REAL(8), INTENT(IN) :: theta_(:)
    integer(4) :: i, j, sid, ierr, icount
    real(8) :: rtemp1, rtemp2, rtemp3, rtemp4, rV1(1)
    integer(4) ::  status(MPI_STATUS_SIZE), iDT_diff(8)

    do i = 1, N_PAR
       IF (theta_(i) .NE. theta_(i)) THEN    ! this will be true if theta_ is NaN
          SMM_distance = NR_BIG2
          return
       endif
    enddo
    if (gcL_partrans) then
       do i = 1, N_PAR
          gr_Vpar(i) = real2ab_exp(gr_paradomain(i,1), gr_paradomain(i,2), theta_(i))
       enddo
    else
       gr_Vpar = theta_
    endif  !if (gcL_partrans)
    
    if (gcI_robust == 21) then
        gr_Vpar(2) = 0.0d0 !no H depreciation at all
    elseif (gcI_robust==12 .OR. gcI_robust==13) then !exo/lbd Hone
        gr_Vpar(3) = 0.0d0 !alpha_I set as 0.0d0
    endif
        
    call sub_node_init()
           
    !=======================================================
    ! broadcast parameters and some state space
    !=======================================================
    if (gcL_debugprint) write(*,'(A,I4,A)') 'iter# ', gcI_itercount, ', Master starts BCAST gr_Vpar'
    
    gr_Vpar_mpi(1:N_PAR) = gr_Vpar
    call MPI_BCAST(gr_Vpar_mpi,N_PAR_mpi,MPI_DOUBLE_PRECISION,master_mpi,MPI_COMM_WORLD,ierr)   
    
    if (gcL_debugprint) write(*,'(A,I4,A,30F20.5)') 'iter# ', gcI_itercount, ', Master done BCAST gr_Vpar=', gr_Vpar
    
    call sub_parse_parameters()
    
    
    !====================================================
    ! Asset space is calculated in the head slaves since
    !   it depends on the heterogenous parameter values.
    !====================================================
    !=======================================================
    ! END of broadcast parameters and some state space
    !=======================================================
    
    !the master is also one of the head slaves
    call cmaxee_VALUEFUNCTION()  !solve DP for given set of parameters (het)
    !start simulation

    if ((gcL_searchinitA) .AND. (gcI_CF==0 .OR. gcI_CF==70 .OR. gcI_CF==88)) then
       pintA_Vmpi(2) = 0.0d0 
       pinitA = max(pinitA_domain(1,1)+1.0d-6, min(pinitA_domain(1,2)*0.5d0, pinitA)) 
       rV1(1) = ab2real_exp(pinitA_domain(1,1), pinitA_domain(1,2), pinitA)
       call Nelder_Meade_C(rV1,1.0d-6,SMM_distanceATinit,1,5000) !Simplex
    else
       rV1(1) = ab2real_exp(pinitA_domain(1,1), pinitA_domain(1,2), pinitA) 
    endif   
    !now the estimated initial asset level
    pintA_Vmpi(2) = 1.0d0
    SMM_distance = SMM_distanceATinit(rV1)
    pinitA = real2ab_exp(pinitA_domain(1,1), pinitA_domain(1,2), rV1(1))

    if (gcL_IsOutput .AND. (gcI_CF==70 .OR. gcI_CF==88)) then
       call DATE_AND_TIME(VALUES=gI_DT1)
       iDT_diff = func_time_elapsed(gI_DT0, gI_DT1)       
       
       write(122,'(A,I4,A,I4,A,I4,A,I4,A,I4,A)') "Master - Time elapsed solving DP & Simulation: ",  &
                 &  iDT_diff(2), " months, ", iDT_diff(3), " days, ", iDT_diff(5), " hrs, ",  &
                 &  iDT_diff(6), " mins, ", iDT_diff(7), " secs."
    endif  !if (gcL_IsOutput)
    
    write(122,'(A,I6,A,F50.10)') "iteration: ", gcI_itercount, ", distance: ", SMM_distance
    
    !!!!! record the minimum distance so far
    if (SMM_distance < grMinDistance) then
       do i = 1, N_PAR
          write(124,*) gr_Vpar(i)
       enddo
       write(124,*) pinitA  !initial asset level
       write(124,'(A,I6,A,F50.10)') "iteration: ", gcI_itercount, ", distance: ", SMM_distance
       grMinDistance = SMM_distance

       ! moments
       if (1>0) then
          if (Ipt==1) then
              i = N_data_0_pt
              j = max(N_data_10, N_data_10_ss, N_data_10_pt)
          else
              i = N_data_0
              j = max(N_data_10, N_data_10_ss)
          endif
          
          do git = i,  j
             if (IHasHealthAge > 0) then
                rtemp1 = MSM_Moments_Sim%a(git,7)  !diff_lfpr_E2G
                rtemp2 = MSM_Moments_Sim%a(git,8)  !diff_lfpr_G2B
                rtemp3 = MSM_Moments_Sim%a(git,9)  !diff_lfpr_B2D
             else
                rtemp1 = -999.0d0
                rtemp2 = -999.0d0
                rtemp3 = -999.0d0
             endif
             
             if (Ipt == 1) then !PT
                 if (IHasHealthAge > 0) then
                     rtemp4 = MSM_Moments_Sim%a(git,10)  !pt
                 else
                     rtemp4 = MSM_Moments_Sim%a(git,7)  !pt
                 endif
             else
                 rtemp4 = -999.0d0
             endif
           
             write(125,'(I8,A,I4,A,F12.4,A,F12.4,A,F12.4,A,F12.4,A,F12.4, &
                   & A,F12.4,A,F12.4,A,F12.4,A,F12.4,A,F12.4,A,F12.4,A,F12.4,A,F12.4)')  &
                   & gcI_itercount, ',', git,         &
                   & ',', MSM_Moments_Sim%a(git,1),   &
                   & ',', MSM_Moments_Sim%a(git,2),   &
                   & ',', MSM_Moments_Sim%a(git,3),   &
                   & ',', MSM_Moments_Sim%a(git,4),   &
                   & ',', MSM_Moments_Sim%a(1,1),     &   !E to U
                   & ',', MSM_Moments_Sim%a(2,1),     &   !U to E
                   & ',', MSM_Moments_Sim%a(git,5),   &   !consumption
                   & ',', MSM_Moments_Sim%a(git,6),   &   !#6 SS
                   & ',', rtemp1, ',', rtemp2, ',', rtemp3, ',', rtemp4,   &   !diff_lfpr_E2G, diff_lfpr_G2B, diff_lfpr_B2D, lfpr_pt
                   & ',', MSM_Moments_Sim%a(3,1)          !wage depreciation after unemployment spell

           enddo
       endif   
    endif
    
    106 continue
    return
  end function SMM_distance
  
  !==========================================================================  
  !distance to be minimized, given the initial Asset level. called in MASTER  
  !==========================================================================
  real(8) function SMM_distanceATinit(Vinit_)
    implicit none
    REAL(8), INTENT(IN) :: Vinit_(:)
    integer(4) :: i, ierr

    do i = 1, 1
       IF (Vinit_(i) .NE. Vinit_(i)) THEN    ! this will be true if Vinit_ is NaN
          SMM_distanceATinit = NR_BIG2
          return
       endif
    enddo
        
    !the initial asset value
    pinitA = real2ab_exp(pinitA_domain(1,1), pinitA_domain(1,2), Vinit_(1))
    !broadcast the initial guess of initial asset to all nodes
    pintA_Vmpi(1) = pinitA
    call MPI_BCAST(pintA_Vmpi,2,MPI_DOUBLE_PRECISION,master_mpi,MPI_COMM_WORLD,ierr)    
    
    !start simulation
    if (gcL_CF_IW_decompose) then
        if (gcI_CF_IW>0) then
            call SMM_SimulateProfiles_IW()
        else
            call SMM_SimulateProfiles()
        endif        
    else
        call SMM_SimulateProfiles()
    endif
    ! Calculate moments
    SMM_distanceATinit = SMM_MomentsDistance()
    
    if (gcL_cal_shadowprice) call SMM_SimulateProfiles_Alt()
    
    return
  end function SMM_distanceATinit
    
  !==============================================================================
  ! Given initial X (H) to Calculate the moments/profiles, and then distance. 
  ! called in all nodes: MASTER, head slaves, and all other slaves.
  !==============================================================================
  real(8) function SMM_MomentsDistance()
    implicit none
    integer(4) :: i, it, j1, j2, j3, j4, j5, ix(3), itemp, tstart, tend,  &
              &   i4n, ierr, iDT_diff(8), jPT, tstart_pt, tend_pt, it_pos
    real(8) :: rw, rh, ev, rtemp, rtemp2, rnullmom
    logical :: lnotpt
    
    
    call MPI_Barrier(MPI_COMM_WORLD,ierr)
    
    tstart_pt = N_data_0_pt  !start age for pt moment
    tend_pt = N_data_10_pt  !end age for pt moment
    if (gcL_IsOutput) then  !periods for simulation: calculate moments for the whole lifecycle
       tstart = N_0
       tend = N_T
    else                    !periods for estimation: calculate moments only for part of lifecycle that we minimize over
       tstart = N_data_0
       tend = max(N_data_10, N_data_10_ss)  !(65,70)
    endif

    !===========================================================================================================================
    ! Calculate moments: 
    ! 1 lfpr_mean, 2 lnw_mean, 3 lnw_sd, 4 lnw_fe_mean, 
    ! 5 consumption  6 ss,  7 diff_lfpr_E2G, 8 diff_lfpr_G2B, 9 diff_lfpr_B2D.
    ! 7 or 10: part-time
    ! 1.1 E to U, 2.1 U to E, 3.1 wage depreciation after unemp
    ! !1-A,2-H,3-AIME,4-c,5-I,6-L,7-Labor,8-Value,9-leisure error
    !===========================================================================================================================
    rnullmom = 0.0d0
    MSM_Moments_Sim%a = 0.0d0
    MSM_sim4m_i%a = 0
    
    ! #7 diff_lfpr_E2G: difference in lfpr for 1 Excellent v.s. 2 Good.
    ! #8 diff_lfpr_G2B: difference in lfpr for 2 Good v.s. 3 Bad.
    ! #9 diff_lfpr_B2D: difference in lfpr for 3 bad v.s. 4 disable.
    if (IHasHealthAge >= 1) then
       j1 = max(N_0+1, N_data_0_health)
       !lfpr for good health and bad health
       do i = gI_Indiv_Tasks_mpi%a(myid_mpi+1,1), gI_Indiv_Tasks_mpi%a(myid_mpi+1,2)
          do it = j1, tend
             it_pos = it - N_0 + 1 
             rtemp = MSM_Indiv%a(6,it_pos,i)  !6-L 
             
             if (Ipt==1 .AND. gcL_skippt .AND. rtemp>0.2d0 .AND. rtemp<0.8d0) cycle !skip working part time
             
             if (MSM_Indiv_Ihealth%a(it_pos,i)==1) then !excellent health 
                MSM_sim4m_i%a(it,1) = MSM_sim4m_i%a(it,1) + 1
                if (rtemp<0.8d0) then
                   MSM_Moments_Sim%a(it,1) = MSM_Moments_Sim%a(it,1) + 1.0d0  !working, excellent health 
                endif
             elseif (MSM_Indiv_Ihealth%a(it_pos,i)==2) then !good health 
                MSM_sim4m_i%a(it,2) = MSM_sim4m_i%a(it,2) + 1
                if (rtemp<0.8d0) then
                   MSM_Moments_Sim%a(it,2) = MSM_Moments_Sim%a(it,2) + 1.0d0  !working, good health
                endif
             elseif (MSM_Indiv_Ihealth%a(it_pos,i)==3) then !bad health 
                MSM_sim4m_i%a(it,3) = MSM_sim4m_i%a(it,3) + 1
                if (rtemp<0.8d0) then
                   MSM_Moments_Sim%a(it,3) = MSM_Moments_Sim%a(it,3) + 1.0d0  !working, bad health
                endif
             elseif (MSM_Indiv_Ihealth%a(it_pos,i)==4) then !disable 
                MSM_sim4m_i%a(it,4) = MSM_sim4m_i%a(it,4) + 1
                if (rtemp<0.8d0) then
                   MSM_Moments_Sim%a(it,4) = MSM_Moments_Sim%a(it,4) + 1.0d0  !working, disable
                endif   
             endif   
          enddo !it
       enddo  !do i
       if (myid_mpi == master_mpi) then  !MASTER
          call MPI_reduce(MPI_In_PLACE,MSM_sim4m_i%a(1,1),N_T*4,MPI_INTEGER,MPI_SUM,master_mpi,MPI_COMM_WORLD,ierr)
          call MPI_reduce(MPI_In_PLACE,MSM_Moments_Sim%a(1,1),N_T*4,MPI_DOUBLE_PRECISION,MPI_SUM,master_mpi,MPI_COMM_WORLD,ierr)
          do it = j1, tend
             if (MSM_sim4m_i%a(it,1) <= 0) then !excellent health
                MSM_Moments_Sim%a(it, 1) = rnullmom  !mean
             else
                MSM_Moments_Sim%a(it, 1) = MSM_Moments_Sim%a(it,1) / dble(MSM_sim4m_i%a(it,1))
             endif
             if (MSM_sim4m_i%a(it,2) <= 0) then !good health
                MSM_Moments_Sim%a(it, 2) = rnullmom  !mean
             else
                MSM_Moments_Sim%a(it, 2) = MSM_Moments_Sim%a(it,2) / dble(MSM_sim4m_i%a(it,2))
             endif
             if (MSM_sim4m_i%a(it,3) <= 0) then !bad health
                MSM_Moments_Sim%a(it, 3) = rnullmom  !mean
             else
                MSM_Moments_Sim%a(it, 3) = MSM_Moments_Sim%a(it,3) / dble(MSM_sim4m_i%a(it,3))
             endif
             if (MSM_sim4m_i%a(it,4) <= 0) then !disable
                MSM_Moments_Sim%a(it, 4) = rnullmom  !mean
             else
                MSM_Moments_Sim%a(it, 4) = MSM_Moments_Sim%a(it,4) / dble(MSM_sim4m_i%a(it,4))
             endif
             MSM_Moments_Sim%a(it, 7) = MSM_Moments_Sim%a(it,1)-MSM_Moments_Sim%a(it,2) !difference: Excellent - Good
             MSM_Moments_Sim%a(it, 8) = MSM_Moments_Sim%a(it,2)-MSM_Moments_Sim%a(it,3) !difference: Good - Bad
             MSM_Moments_Sim%a(it, 9) = MSM_Moments_Sim%a(it,3)-MSM_Moments_Sim%a(it,4) !difference: Bad - Disable
          enddo  !do it = tstart, tend
       else
          call MPI_reduce(MSM_sim4m_i%a(1,1),MSM_sim4m_i%a(1,1),N_T*4,MPI_INTEGER,MPI_SUM,master_mpi,MPI_COMM_WORLD,ierr)
          call MPI_reduce(MSM_Moments_Sim%a(1,1),MSM_Moments_Sim%a(1,1),N_T*4, &
                       &  MPI_DOUBLE_PRECISION,MPI_SUM,master_mpi,MPI_COMM_WORLD,ierr) 
       endif  !if (myid_mpi == master_mpi)
       MSM_Moments_Sim%a(:,1:4) = 0.0d0
       MSM_sim4m_i%a(:,1:4) = 0
    endif !if (IHasHealthAge >= 1) then    
    
    !=== 7 or 10 part time
    if (Ipt == 1) then
       if (IHasHealthAge >= 1) then 
           jPT = 10
       else
           jPT = 7
       endif
       
       do i = gI_Indiv_Tasks_mpi%a(myid_mpi+1,1), gI_Indiv_Tasks_mpi%a(myid_mpi+1,2)
          do it = tstart_pt, tend_pt   !tstart, tend
             it_pos = it - N_0 + 1 
             !=== 7/10 LFPR - part time
             MSM_sim4m_i%a(it,jPT) = MSM_sim4m_i%a(it,jPT) + 1
             if (MSM_Indiv%a(6,it_pos,i)>=0.2d0 .AND. MSM_Indiv%a(6,it_pos,i)<=0.8d0) then
                MSM_Moments_Sim%a(it,jPT) = MSM_Moments_Sim%a(it,jPT) + 1.0d0 !working part time
             endif
          enddo !it
       enddo
       if (myid_mpi == master_mpi) then  !MASTER
          call MPI_reduce(MPI_In_PLACE,MSM_sim4m_i%a(1,jPT),N_T,MPI_INTEGER,MPI_SUM,master_mpi,MPI_COMM_WORLD,ierr)
          call MPI_reduce(MPI_In_PLACE,MSM_Moments_Sim%a(1,jPT),N_T,MPI_DOUBLE_PRECISION,MPI_SUM,master_mpi,MPI_COMM_WORLD,ierr)
          do it = tstart_pt, tend_pt
             if (MSM_sim4m_i%a(it,jPT) <= 0) then
                MSM_Moments_Sim%a(it, jPT) = rnullmom  !mean
             else
                MSM_Moments_Sim%a(it, jPT) = MSM_Moments_Sim%a(it,jPT) / dble(MSM_sim4m_i%a(it,jPT))
             endif
          enddo  !do it = tstart, tend
       else
          call MPI_reduce(MSM_sim4m_i%a(1,jPT),MSM_sim4m_i%a(1,jPT),N_T,MPI_INTEGER,MPI_SUM,master_mpi,MPI_COMM_WORLD,ierr)
          call MPI_reduce(MSM_Moments_Sim%a(1,jPT),MSM_Moments_Sim%a(1,jPT),N_T, &
                       &  MPI_DOUBLE_PRECISION,MPI_SUM,master_mpi,MPI_COMM_WORLD,ierr) 
       endif  !if (myid_mpi == master_mpi)
    endif !if (Ipt == 1) then 
    

    !in each slaves (head and others)
    do i = gI_Indiv_Tasks_mpi%a(myid_mpi+1,1), gI_Indiv_Tasks_mpi%a(myid_mpi+1,2)
       do it = tstart, tend
          it_pos = it - N_0 + 1 
          
          !=== 1 labor force participation rate 
          if (Ipt==1 .AND. gcL_skippt .AND. MSM_Indiv%a(6,it_pos,i)>0.2d0 .AND. MSM_Indiv%a(6,it_pos,i)<0.8d0) then  !skip working part time
              lnotpt = .FALSE.
          else
              lnotpt = .TRUE.
          endif
          
          if (lnotpt) MSM_sim4m_i%a(it,1) = MSM_sim4m_i%a(it,1) + 1     !1 lfpr: no part time
          
          if (MSM_Indiv%a(6,it_pos,i)<0.9d0) then   !working, full time or part time
              if (lnotpt) MSM_Moments_Sim%a(it,1) = MSM_Moments_Sim%a(it,1) + 1.0d0  !1 lfpr: no part time
              
              !=== 2 lnw_mean (incl pt), 4 lnw_fd_mean (excl pt, was lnw_fe_mean)
              rtemp2 = pw(it) * MSM_Indiv%a(2,it_pos,i)*MSM_Indiv%a(7,it_pos,i)      !actual labor income
              rtemp = 1.0d0 - MSM_Indiv%a(6,it_pos,i)    !observed labor. always 1 in discrete labor case.
              rtemp = log(rtemp2 / rtemp)  !log wage
              if (rtemp>grMinlnw )  then   !actual lnw > grMinlnw
                  MSM_sim4m_i%a(it,2) = MSM_sim4m_i%a(it,2) + 1               !2 lnw_mean, including PT
                  MSM_Moments_Sim%a(it,2) = MSM_Moments_Sim%a(it,2) + rtemp   !2 lnw_mean, including PT
             
                  !for #4 fixed effects, from N_data_0 to N_data_10, (FD) from FT to FT only.
                  if (lnotpt .AND. (it>=MSM_Indiv_Ihet%a(i,2)) .and. (it<=MSM_Indiv_Ihet%a(i,3))) then  !working FT at t and observed
                      if (it>N_data_0 .AND. (MSM_Indiv%a(6,it_pos-1,i)<0.1d0)) then  !working FT at t-1
                          rtemp2 = pw(it-1) * MSM_Indiv%a(2,it_pos-1,i)*MSM_Indiv%a(7,it_pos-1,i)   !actual labor income at t-1
                          rtemp2 = log( rtemp2 / (1.0d0-MSM_Indiv%a(6,it_pos-1,i)) )  !lnw at t-1
                          if (rtemp2>grMinlnw) then
                              MSM_sim4m_i%a(it,4) = MSM_sim4m_i%a(it,4) + 1                      !4 lnw_fd
                              MSM_Moments_Sim%a(it,4) = MSM_Moments_Sim%a(it,4) + rtemp-rtemp2   !dlnw
                              if (it==N_data_0+1) then
                                  MSM_sim4m_i%a(it-1,4) = MSM_sim4m_i%a(it-1,4) + 1                !lnw at t-1, serving as lnw_fd_base
                                  MSM_Moments_Sim%a(it-1,4) = MSM_Moments_Sim%a(it-1,4) + rtemp2   !lnw at t-1
                              endif
                          endif  !if (rtemp2>grMinlnw)
                      endif   !if (it>N_data_0)
                  endif  !  if (it<=N_data_10) 
             
                  !for 3.1 wage depreciation after unemp, between 41 and 65 (both inclusive)
                  if (((it>40) .and. (it<=N_data_10)) .AND.  &   !in data we use the +2 waves (3 months) as the reference
                      &  ((it>=MSM_Indiv_Ihet%a(i,2)) .and. (it<=MSM_Indiv_Ihet%a(i,3))) ) then
                       if ((MSM_Indiv%a(6,it_pos-2,i)<0.1d0) .AND. (MSM_Indiv%a(6,it_pos-1,i)>0.9d0) .AND.  &
                         & (MSM_Indiv%a(6,it_pos,i)<0.1d0)) then  !it-2: working FT, it-1 not working, it: working FT
                           rtemp2 = pw(it-2) * MSM_Indiv%a(2,it_pos-2,i)*MSM_Indiv%a(7,it_pos-2,i)   !actual labor income at t-2
                           rtemp2 = log( rtemp2 / (1.0d0-MSM_Indiv%a(6,it_pos-2,i)) )  !lnw at t-2
                           if (rtemp2>grMinlnw) then
                              MSM_sim4m_i%a(3,1) = MSM_sim4m_i%a(3,1) + 1                          
                              MSM_Moments_Sim%a(3,1) = MSM_Moments_Sim%a(3,1) + rtemp - rtemp2
                           endif  !if (rtemp2>grMinlnw)
                       endif     !if ((MSM_Indiv%a(6,it_pos-1,i)>0.9d0)
                  endif    !if (it>40) 
              endif   !if ( rtemp>grMinlnw )
          endif  !if (MSM_Indiv%a(6,it_pos,i)<0.9d0)
          
                    
          !=== 5 consumption
          MSM_sim4m_i%a(it,5) = MSM_sim4m_i%a(it,5) + 1
          if (MSM_Indiv_IFS%a(it_pos,i) == 1) then  !Family status
             j3 = 1
          else
             j3 = 2 
          endif
          MSM_Moments_Sim%a(it,5) = MSM_Moments_Sim%a(it,5)  &
               & + MSM_Indiv%a(4,it_pos,i)/grAdultEquivScale(j3)  !adult equivalent consumption
          
          ! ===== #6. SS (ssa), from N_data_0_ss (62) to N_data_10_ss (69)
          !MSM_Indiv_Irecss%a(it_pos,i) is the status coming to it
          !MSM_Indiv_Irecss%a(it_pos+1,i) is the decision at it. 1 means receiving SS at it, NOT it+1
          MSM_sim4m_i%a(it,6) = MSM_sim4m_i%a(it,6) + 1
          if ((it<N_T) .and. (MSM_Indiv_Irecss%a(it_pos,i)==0)) then
             if (MSM_Indiv_Irecss%a(it_pos+1,i)==1) MSM_Moments_Sim%a(it,6)=MSM_Moments_Sim%a(it,6)+1.0d0
          endif
          
          !=== 1.1 E to U, 2.1 U to E. Aggregate between 35 and 50 (both inclusive). SKIP PT
          if (lnotpt .AND. (it>=35 .AND. it<=50) .AND.  &   !aggregate over 35-50
            &  ((it>=MSM_Indiv_Ihet%a(i,2)) .and. (it+1<=MSM_Indiv_Ihet%a(i,3))) ) then              
              if (MSM_Indiv%a(6,it_pos+1,i)>0.2d0 .AND. MSM_Indiv%a(6,it_pos+1,i)<0.8d0) then  !t+1: In SIPP we skip PT.
              else
                  if (MSM_Indiv%a(6,it_pos,i) < 0.9d0) then  !working (FT), E now, at each age
                     MSM_sim4m_i%a(1,1) = MSM_sim4m_i%a(1,1) + 1
                     if (MSM_Indiv%a(6,it_pos+1,i) > 0.9d0) then
                        MSM_Moments_Sim%a(1,1) = MSM_Moments_Sim%a(1,1) + 1.0d0   !E to U, aggregate
                     endif
                  else  !not working, U now
                     MSM_sim4m_i%a(2,1) = MSM_sim4m_i%a(2,1) + 1
                     if (MSM_Indiv%a(6,it_pos+1,i) < 0.9d0) then
                        MSM_Moments_Sim%a(2,1) = MSM_Moments_Sim%a(2,1) + 1.0d0   !U to E, aggregate
                     endif
                  endif !if (MSM_Indiv%a(6,it_pos,i) < 0.9d0)
              endif  !if (Ipt==1 .AND. gcL_skippt .AND.    
          endif  !if (it>=35 .AND. it<=50)
          
       enddo !it
    enddo  !do i
    !=== #1 labor force participation rate, and transitions (1.1., 2.1.)
    !=== #4 lnw_fd, in the master only
    !=== #5 consumption, in the master only
    !=== #6. frSS (ssa)
    if (myid_mpi == master_mpi) then
       call MPI_Reduce(MPI_In_PLACE,MSM_sim4m_i%a(1,1),N_T,MPI_INTEGER,MPI_SUM,master_mpi,MPI_COMM_WORLD,ierr)
       call MPI_reduce(MPI_In_PLACE,MSM_Moments_Sim%a(1,1),N_T,MPI_DOUBLE_PRECISION,MPI_SUM,master_mpi,MPI_COMM_WORLD,ierr)
       call MPI_Reduce(MPI_In_PLACE,MSM_sim4m_i%a(1,4),N_T*3,MPI_INTEGER,MPI_SUM,master_mpi,MPI_COMM_WORLD,ierr)
       call MPI_reduce(MPI_In_PLACE,MSM_Moments_Sim%a(1,4),N_T*3,MPI_DOUBLE_PRECISION,MPI_SUM,master_mpi,MPI_COMM_WORLD,ierr)
    else
       call MPI_Reduce(MSM_sim4m_i%a(1,1),MSM_sim4m_i%a(1,1),N_T,MPI_INTEGER,MPI_SUM,master_mpi,MPI_COMM_WORLD,ierr) 
       call MPI_reduce(MSM_Moments_Sim%a(1,1),MSM_Moments_Sim%a(1,1),N_T,MPI_DOUBLE_PRECISION,MPI_SUM, &
                    &  master_mpi,MPI_COMM_WORLD,ierr) 
       call MPI_Reduce(MSM_sim4m_i%a(1,4),MSM_sim4m_i%a(1,4),N_T*3,MPI_INTEGER,MPI_SUM,master_mpi,MPI_COMM_WORLD,ierr) 
       call MPI_reduce(MSM_Moments_Sim%a(1,4),MSM_Moments_Sim%a(1,4),N_T*3,MPI_DOUBLE_PRECISION,MPI_SUM, &
                    &  master_mpi,MPI_COMM_WORLD,ierr) 
    endif
            
    if (myid_mpi == master_mpi) then !MASTER
       do it = tstart, tend
          if (MSM_sim4m_i%a(it,1) <= 0) then !#1 lfpr
             MSM_Moments_Sim%a(it, 1) = rnullmom  !mean
          else
             MSM_Moments_Sim%a(it, 1) = MSM_Moments_Sim%a(it,1) / dble(MSM_sim4m_i%a(it,1))
          endif
          if (MSM_sim4m_i%a(it,5) <= 0) then  !#5 consumption
             MSM_Moments_Sim%a(it, 5) = rnullmom  !mean
          else
             MSM_Moments_Sim%a(it, 5) = MSM_Moments_Sim%a(it,5) / dble(MSM_sim4m_i%a(it,5))
          endif
          if (MSM_sim4m_i%a(it,6) <= 0) then  !#6 SS
             MSM_Moments_Sim%a(it, 6) = rnullmom  !mean
          else
             MSM_Moments_Sim%a(it, 6) = MSM_Moments_Sim%a(it,6) / dble(MSM_sim4m_i%a(it,6))
          endif
       enddo  !do it = tstart, tend
       
       do i = 1, 1  !aggregate transition 1.1, 2.1, and wage depreciation 3.1
          do it = 1, 3
             if (MSM_sim4m_i%a(it,i) <= 0) then
                MSM_Moments_Sim%a(it, i) = rnullmom  !mean
             else
                MSM_Moments_Sim%a(it, i) = MSM_Moments_Sim%a(it,i) / dble(MSM_sim4m_i%a(it,i))
             endif
          enddo !do it
       enddo  !do i = 1,1
    endif  !if (myid_mpi == master_mpi)
    
    !=== 2 lnw_mean, in all nodes, needed in all nodes for lnw_sd
    call MPI_Allreduce(MPI_In_PLACE,MSM_sim4m_i%a(1,2),N_T,MPI_INTEGER,MPI_SUM,MPI_COMM_WORLD,ierr)
    call MPI_Allreduce(MPI_In_PLACE,MSM_Moments_Sim%a(1,2),N_T,MPI_DOUBLE_PRECISION,MPI_SUM,MPI_COMM_WORLD,ierr)
    do it = tstart, tend
       if (MSM_sim4m_i%a(it,2) <= 0) then
          MSM_Moments_Sim%a(it, 2) = rnullmom  !mean
       else
          MSM_Moments_Sim%a(it, 2) = MSM_Moments_Sim%a(it,2) / dble(MSM_sim4m_i%a(it,2))
       endif
    enddo  !do it = tstart, tend
    
    !#4 lnw_fd_mean
    if (myid_mpi == master_mpi) then !MASTER
       do it = N_data_0, N_data_10
          if (MSM_sim4m_i%a(it,4) <= 0) then 
             MSM_Moments_Sim%a(it, 4) = rnullmom  !mean
          else
             MSM_Moments_Sim%a(it, 4) = MSM_Moments_Sim%a(it,4) / dble(MSM_sim4m_i%a(it,4))
          endif
          if (it>N_data_0) MSM_Moments_Sim%a(it, 4) = MSM_Moments_Sim%a(it, 4) + MSM_Moments_Sim%a(it-1, 4)
       enddo !do it = N_data_0, N_data_10
    endif  !if (myid_mpi == master_mpi)
    
    !=== 3 lnw_sd
    do i = gI_Indiv_Tasks_mpi%a(myid_mpi+1,1), gI_Indiv_Tasks_mpi%a(myid_mpi+1,2)
       do it = tstart, tend
          it_pos = it - N_0 + 1 
          if (MSM_Indiv%a(6,it_pos,i)<0.9d0)  then   !leisure
             rtemp = 1.0d0 - MSM_Indiv%a(6,it_pos,i)   !observed labor. always 1 in discrete labor case.
             rtemp2 = pw(it) * MSM_Indiv%a(2,it_pos,i)*MSM_Indiv%a(7,it_pos,i)            !actual labor income
             rtemp = log(rtemp2 / rtemp)  !log wage
             if (rtemp>grMinlnw) then
                 MSM_Moments_Sim%a(it,3) = MSM_Moments_Sim%a(it,3) + (rtemp-MSM_Moments_Sim%a(it,2))**2 !log wage
             endif
          endif  !if ( (rtemp2>grMinlnw) .and. (MSM_Indiv%a(6,it,i)<1.0d0))
       enddo !it
    enddo  !do i
    if (myid_mpi == master_mpi) then
       call MPI_reduce(MPI_In_PLACE,MSM_Moments_Sim%a(1,3),N_T,MPI_DOUBLE_PRECISION, &
                    &  MPI_SUM,master_mpi,MPI_COMM_WORLD,ierr)
    else
       call MPI_reduce(MSM_Moments_Sim%a(1,3),MSM_Moments_Sim%a(1,3),N_T,MPI_DOUBLE_PRECISION, &
                    &  MPI_SUM,master_mpi,MPI_COMM_WORLD,ierr) 
    endif
    if (myid_mpi == master_mpi) then !MASTER
       do it = tstart, tend
          if (MSM_sim4m_i%a(it,2) <= 1) then
             MSM_Moments_Sim%a(it, 3) = rnullmom  !mean
          else
             MSM_Moments_Sim%a(it, 3) = sqrt(MSM_Moments_Sim%a(it,3) / dble(MSM_sim4m_i%a(it,2)-1))
          endif
       enddo  !do it = tstart, tend
    endif  !if (myid_mpi == master_mpi)
            
        
    !===========================================================================================================================
    ! Calculate moments: 
    ! 1 lfpr_mean, 2 lnw_mean, 3 lnw_sd, 4 lnw_fd_mean, 
    ! 5 consumption  6 ss,  7 diff_lfpr_E2G, 8 diff_lfpr_G2B, 9 diff_lfpr_B2D
    ! 7 or 10: part-time
    ! 1.1 E to U, 2.1 U to E, 3.1 wage depreciation after unemp
    !===========================================================================================================================
    SMM_MomentsDistance = 0.0d0
    i = 0 !second dimension, column i of the weight matrix MSM_wgt%a(ix(1), i)
    if (myid_mpi == master_mpi) then
       do j2 = 1, N_moments_type
          j4 = N_data_10 
          if (j2 == 5) then  !5 consumption
             j3 = N_data_0_c
          elseif (j2 == 6) then  !6 ss
             j3 = N_data_0_ss
             j4 = N_data_10_ss
          elseif ((j2>=7) .AND. (j2<=10)) then  !#7 #8 #9 diff_lfpr and/or 7/10 PT
             if (IHasHealthAge>0 .AND. ((j2>=7) .AND. (j2<=9)) ) then
                j3 = N_data_0_health       !diff_lfpr
             elseif (Ipt==1 .and. (j2==7 .OR. j2==10)) then  !pt
                j3 = N_data_0_pt
                j4 = N_data_10_pt
             else
                cycle 
             endif
          else
             j3 = N_data_0 
          endif
          do j1 = j3, j4
             !column i (mapping (j1,j2) in MSM_Moments_Sim) of the weight matrix MSM_wgt%a(ix(1), i)
             i = i + 1
             !!!! first dimension now
             rtemp = 0.0d0
             ix(1) = 0  !first dimension
             do itemp = 1, N_moments_type
                j5 = N_data_10 
                if (itemp == 5) then  !#5 consumption
                   ix(3) = N_data_0_c
                elseif (itemp == 6) then  !#6 SS
                   ix(3) = N_data_0_ss
                   j5 = N_data_10_ss
                elseif ((itemp>=7) .AND. (itemp<=10)) then  !#7 #8 #9 diff_lfpr and/or 7/10 PT
                   if (IHasHealthAge>0 .AND. ((itemp>=7) .AND. (itemp<=9))) then
                      ix(3) = N_data_0_health  !diff_lfpr
                   elseif (Ipt==1 .and. (itemp==7 .OR. itemp==10)) then
                      ix(3) = N_data_0_pt  !pt
                      j5 = N_data_10_pt
                   else
                      cycle
                   endif
                else
                   ix(3) = N_data_0 
                endif  !if (itemp == N_moments_type)
                do it = ix(3), j5
                   ix(1) = ix(1) + 1 
                   if (MSM_Moments_Data%a(it,itemp) > -900.0d0)  &
                     &  rtemp = rtemp + (MSM_Moments_Sim%a(it,itemp)-MSM_Moments_Data%a(it,itemp)) * MSM_wgt%a(ix(1), i)
                enddo !itemp
             enddo !itemp
             !!!! end of first dimension
             if (MSM_Moments_Data%a(j1,j2) > -900.0d0) then 
                rtemp2 = rtemp * (MSM_Moments_Sim%a(j1,j2)-MSM_Moments_Data%a(j1,j2))
             else
                rtemp2 = 0.0d0
             endif
             SMM_MomentsDistance = SMM_MomentsDistance + rtemp2
          enddo !j1          
       enddo ! j2
       !1.1 E to U, 2.1 U to E, 3.1 lnw depreciation
       do j2 = 1, 3
          j1 = N_moments_all - 3 + j2 
          SMM_MomentsDistance = SMM_MomentsDistance & 
               &  + ((MSM_Moments_Sim%a(j2,1)-MSM_Moments_Data%a(j2,1))**2) * MSM_wgt%a(j1, j1)
       enddo !do j2
       SMM_MomentsDistance = SMM_MomentsDistance * dble(N_Data) / (1.0d0+dble(N_Data)/dble(N_Sim))  !Chi^2
    endif !if (myid_mpi == master_mpi)
        
    if (gcL_IsOutput) then  !to print out simulated individual results
       i = min(gI_Indiv_Tasks_mpi%a(myid_mpi+1,1), N_Sim) 
       call MPI_Gatherv(MSM_Indiv%a(1,1,i), gI_Indiv_size_mpi*gI_Indiv_Tasks_mpi%a(myid_mpi+1,3), MPI_DOUBLE_PRECISION,   & 
                     &  MSM_Indiv%a(1,1,1), gI_Indiv_size_mpi*gI_Indiv_Tasks_mpi%a(:,3), &
                     &  gI_Indiv_size_mpi*gI_Indiv_Tasks_mpi%a(:,4),     &
                     &  MPI_DOUBLE_PRECISION, master_mpi, MPI_COMM_WORLD, ierr)
       call MPI_Gatherv(MSM_Indiv_Irecss%a(1,i), N_T_grid*gI_Indiv_Tasks_mpi%a(myid_mpi+1,3), MPI_INTEGER,   & 
                     &  MSM_Indiv_Irecss%a(1,1), N_T_grid*gI_Indiv_Tasks_mpi%a(:,3),  &
                     &  N_T_grid*gI_Indiv_Tasks_mpi%a(:,4),     &
                     &  MPI_INTEGER, master_mpi, MPI_COMM_WORLD, ierr)
    endif  !if (gcL_IsOutput)
    
    return
    
  end function SMM_MomentsDistance
  
  
  !========================================================================================================
  ! Simulate lifecycle profiles all by interpolation, called by all nodes (master, head slaves and slaves)
  !========================================================================================================
  subroutine SMM_SimulateProfiles()
    implicit none
    integer(4) :: i, j, ix(3), ierr
    real(8) :: rtemp, rtemp2, rcx3(3), rcx2(2), ev, rVy(3), rptx
    
    !MSM_Indiv
    !1-A,2-H,3-AIME,4-c,5-I,6-L,7-Labor,8-Value,9-leisure error,10-XI error,
    !11-HEALTH shock,12-V0,13-V1, 14-Prob of leisure=0 or Phi(epsilon_tmin) in PT, 
    !15-FS error, 16-FS spouse inc
    !17-produced output, 18-ss tax, 19-medicare tax, 20-income tax, 21-received ss benefit, 22-ssdiinc
    !23-V(-1), 24-Phi(epsilon_tmax)
    
    if (gcI_CF == -12) then      !-12 ---Elasticity, unanticipated
        pw(gcI_elas_t) = 1.1d0
    endif  !if (gcI_CF == -12)    
    
    !We need to start from N_0 where AIME=0
    MSM_Indiv%a(8,:,:) = - NR_BIG2 !VALUE
    MSM_Indiv%a(17:22,:,:) = 0.0d0  !17-produced output, 18-ss tax, 19-medicare tax, 20-income tax, 21-received ss benefit, 22-ssdiinc
    MSM_Indiv_Irecss%a(:,:) = 0    !Reset RECSS
    do i = gI_Indiv_Tasks_mpi%a(myid_mpi+1,1), gI_Indiv_Tasks_mpi%a(myid_mpi+1,2) 
       !========= Initial conditions ===============
       git = N_0
       git_pos = git - N_0 + 1
       MSM_Indiv%a(1,git_pos,i) = pinitA + MSM_Indiv_initX%a(i,1)      !initial Asset
       !MSM_Indiv%a(1,git_pos,i) = max(sgridA%a(git_pos,1), min(sgridA%a(git_pos,VN_A(git)), MSM_Indiv%a(1,git_pos,i))) !initial A
          
       MSM_Indiv%a(2,git_pos,i) = pinitHmean + pinitHsd * MSM_Indiv_initX%a(i,2) !initial H
       
       if (N_het_H0>1) then
          MSM_Indiv%a(2,git_pos,i) = MSM_Indiv%a(2,git_pos,i) + pinitH_a0*pa0 + pinitH_pi*ppi
       endif
       
       MSM_Indiv%a(2,git_pos,i) = exp(MSM_Indiv%a(2,git_pos,i))   !initial H
       
       MSM_Indiv%a(2,git_pos,i) = max(sgridH%a(git_pos,1)+1.0d-6, min(sgridH%a(git_pos,N_H)-1.0d-6, MSM_Indiv%a(2,git_pos,i))) !H
       !MSM_Indiv%a(3,git_pos,i) = sgridAIME%a(git_pos,INT(sgridAIME%d2/2)) !AIME
       MSM_Indiv%a(3,git_pos,i) = min(sgridAIME%a(git_pos,1)+1.0d-6, sgridAIME%a(git_pos,VN_AIME(git)))  !initial AIME at N_0
       
       !1 RECSS 2 HEALTH
       MSM_Indiv_Irecss%a(git_pos,i) = 0  !this is the status of current period: not receive ss at N_0
       
       !health status at N_0 is set at SUBROUTINE INIT_HEALTH
       !END of ========= Initial conditions ===============
       
       do git = N_0, N_T
          git_pos = git - N_0 + 1 
          call xPosadjust(git, MSM_Indiv%a(1:3,git_pos,i), ix)
          gA = MSM_Indiv%a(1,git_pos,i)
          gH = MSM_Indiv%a(2,git_pos,i)
          gAIME = MSM_Indiv%a(3,git_pos,i)
          giA = ix(1)
          giH = ix(2)
          giAIME = ix(3)
          giRECSS = MSM_Indiv_Irecss%a(git_pos,i) + 1   !index, 1-not RECSS, 2-RECSS
          if (giRECSS<1 .OR. giRECSS>2) then
             write(*,*) 'WRONG 1: i=',i,',git=',git,',giRECSS=',giRECSS, &
                      & ',MSM_Indiv_Irecss%a(N_0:N_T,i)=',MSM_Indiv_Irecss%a(N_0-N_0+1:N_T-N_0+1,i)
          endif
          giHEALTH = MSM_Indiv_Ihealth%a(git_pos,i)      !health status, 1-excellent, 2-good, 3-bad 4-disabled
          giFS = MSM_Indiv_IFS%a(git_pos,i)   !Family status: 1-3
          
          if (gcI_CF == -12 .AND. git == gcI_elas_t) then  !!unanticipated elasticity calculation. Use optimization when shock happens
              call SMM_SimulateProfiles_optimization(i)
          else   !use interpolation
              !First, get depsilon to decide leisure being 0 or 1
              if (Ipt == 0)  then
                 if (gcL_triangle) then
                    rtemp = linear_inte_tetrahedra3(sgridA%a(git_pos,ix(1):(ix(1)+1)), &
                         &  sgridAIME%a(git_pos,ix(3):(ix(3)+1)), sgridH%a(git_pos,ix(2):(ix(2)+1)), &
                         &  depsilon(0)%a(ix(1):(ix(1)+1),ix(3):(ix(3)+1),giHEALTH,giRECSS,giFS,ix(2):(ix(2)+1),git_pos),  &
                         &  (/gA,gAIME,gH/) ) !threshold
                 else
                    rtemp = linear_inte3(gcI_inte_IL_log, (/-1,-1,1/), sgridA%a(git_pos,ix(1):(ix(1)+1)), &
                         &  sgridAIME%a(git_pos,ix(3):(ix(3)+1)), sgridH%a(git_pos,ix(2):(ix(2)+1)), &
                         &  depsilon(0)%a(ix(1):(ix(1)+1),ix(3):(ix(3)+1),giHEALTH,giRECSS,giFS,ix(2):(ix(2)+1),git_pos),  &
                         &  (/gA,gAIME,gH/) ) !threshold
                 endif
              else   !(Ipt == 1)
                                    
                  !!!! calculate V0, V1, and Vp, and then derive two thresholds.
                  !!assuming there is always stochastic shock
                  do j = 0, Ipt+1
                     rVy(1+j) = linear_inte_tetrahedra3(sgridA%a(git_pos,ix(1):(ix(1)+1)), &
                         &  sgridAIME%a(git_pos,ix(3):(ix(3)+1)), sgridH%a(git_pos,ix(2):(ix(2)+1)), &
                         &  dvV(j)%a(ix(1):(ix(1)+1),ix(3):(ix(3)+1),giHEALTH,giRECSS,giFS,ix(2):(ix(2)+1),git_pos),  &
                         &  (/gA,gAIME,gH/) ) !v0
                  enddo !do j
                  rcx3(1) = rVy(1)  !leisure==0
                  rcx3(2) = rVy(2)  !leisure==1
                  rptx = rVy(3)     !leisure==2 (pt)
                  
                  rVy(3) = gr_leisure_taste(git,giHEALTH,giFS)
                  rcx3(3) = rcx3(1) - rcx3(2)  !leisure 0 - leisure 1
                  rVy(1) = min((rcx3(1)-rptx)/max(NR_EPS,plpt_t(git,giFS)), rcx3(3))   !gamma_tmin. Condition A8. Leisure 0 - p
                  rVy(2) = max((rptx-rcx3(2))/max(NR_EPS,(1.0d0-plpt_t(git,giFS))), rcx3(3))  !gamma_tmax. Leisure p - 1
                                     
                  rtemp2 = (log(max(NR_EPS, rVy(1)))-rVy(3)-grUVmultiplier_ln)/pa_epsilon !epsilon_tmin. /pa_epsilon is to get a standard normal with sd=1
                  rtemp = (log(max(NR_EPS, rVy(2)))-rVy(3)-grUVmultiplier_ln)/pa_epsilon !epsilon_tmax
              endif  !if (Ipt == 0)
                            
              if (Ipt == 0)  then  !full time only
                  if (MSM_Indiv%a(9,git_pos,i) <= rtemp) then !leisure = 0
                      MSM_Indiv%a(6,git_pos,i) = 0.0d0
                  else
                      MSM_Indiv%a(6,git_pos,i) = 1.0d0
                  endif
                  MSM_Indiv%a(14,git_pos,i) = CDF_Normal(rtemp, 0.0d0, 1.0d0) !prob of leisure=0
              else   !Ipt == 1, part time
                  if (MSM_Indiv%a(9,git_pos,i) <= rtemp2) then !leisure = 0
                      MSM_Indiv%a(6,git_pos,i) = 0.0d0
                  elseif (MSM_Indiv%a(9,git_pos,i) <= rtemp) then !leisure = 0.5
                      MSM_Indiv%a(6,git_pos,i) = 0.5d0
                  else
                      MSM_Indiv%a(6,git_pos,i) = 1.0d0
                  endif
                  MSM_Indiv%a(14,git_pos,i) = CDF_Normal(rtemp2, 0.0d0, 1.0d0)  !prob of leisure=0
                  MSM_Indiv%a(24,git_pos,i) = CDF_Normal(rtemp, 0.0d0, 1.0d0)   !prob of leisure=pt or 0
              endif  !if (Ipt == 0)
              
              if (gcL_triangle) then          
                 MSM_Indiv%a(12,git_pos,i) = linear_inte_tetrahedra3(sgridA%a(git_pos,ix(1):(ix(1)+1)), &
                         &  sgridAIME%a(git_pos,ix(3):(ix(3)+1)), sgridH%a(git_pos,ix(2):(ix(2)+1)), &
                         &  dvV(0)%a(ix(1):(ix(1)+1),ix(3):(ix(3)+1),giHEALTH,giRECSS,giFS,ix(2):(ix(2)+1),git_pos),  &
                         &  (/gA,gAIME,gH/) ) !V0
                 MSM_Indiv%a(13,git_pos,i) = linear_inte_tetrahedra3(sgridA%a(git_pos,ix(1):(ix(1)+1)), &
                         &  sgridAIME%a(git_pos,ix(3):(ix(3)+1)), sgridH%a(git_pos,ix(2):(ix(2)+1)), &
                         &  dvV(1)%a(ix(1):(ix(1)+1),ix(3):(ix(3)+1),giHEALTH,giRECSS,giFS,ix(2):(ix(2)+1),git_pos),  &
                         &  (/gA,gAIME,gH/) ) !V1
                 if (Ipt==1) then
                     MSM_Indiv%a(23,git_pos,i) = linear_inte_tetrahedra3(sgridA%a(git_pos,ix(1):(ix(1)+1)), &
                         &  sgridAIME%a(git_pos,ix(3):(ix(3)+1)), sgridH%a(git_pos,ix(2):(ix(2)+1)), &
                         &  dvV(2)%a(ix(1):(ix(1)+1),ix(3):(ix(3)+1),giHEALTH,giRECSS,giFS,ix(2):(ix(2)+1),git_pos),  &
                         &  (/gA,gAIME,gH/) ) !V(2)
                 endif
              else
                 MSM_Indiv%a(12,git_pos,i) = linear_inte3(gcI_inte_v_log, (/1,1,1/), sgridA%a(git_pos,ix(1):(ix(1)+1)), &
                         &  sgridAIME%a(git_pos,ix(3):(ix(3)+1)), sgridH%a(git_pos,ix(2):(ix(2)+1)), &
                         &  dvV(0)%a(ix(1):(ix(1)+1),ix(3):(ix(3)+1),giHEALTH,giRECSS,giFS,ix(2):(ix(2)+1),git_pos),  &
                         &  (/gA,gAIME,gH/) ) !V0
                 MSM_Indiv%a(13,git_pos,i) = linear_inte3(gcI_inte_v_log, (/1,1,1/), sgridA%a(git_pos,ix(1):(ix(1)+1)), &
                         &  sgridAIME%a(git_pos,ix(3):(ix(3)+1)), sgridH%a(git_pos,ix(2):(ix(2)+1)), &
                         &  dvV(1)%a(ix(1):(ix(1)+1),ix(3):(ix(3)+1),giHEALTH,giRECSS,giFS,ix(2):(ix(2)+1),git_pos),  &
                         &  (/gA,gAIME,gH/) ) !V1
                 if (Ipt==1) then
                     MSM_Indiv%a(23,git_pos,i) = linear_inte3(gcI_inte_v_log, (/1,1,1/), sgridA%a(git_pos,ix(1):(ix(1)+1)), &
                         &  sgridAIME%a(git_pos,ix(3):(ix(3)+1)), sgridH%a(git_pos,ix(2):(ix(2)+1)), &
                         &  dvV(2)%a(ix(1):(ix(1)+1),ix(3):(ix(3)+1),giHEALTH,giRECSS,giFS,ix(2):(ix(2)+1),git_pos),  &
                         &  (/gA,gAIME,gH/) ) !V(2)
                 endif
              endif 
              
              if (giRECSS == 2) then  !already receiving ss 
                 !the filter makes sure that now git >= pERA
                 if (git < N_T) MSM_Indiv_Irecss%a(git_pos+1,i) = 1  !receiving ss is absorbing
                 giRECSS_next = 2
                 call sub_ForwardInt(MSM_Indiv%a(4:6,git_pos,i),rcx2,MSM_Indiv%a(8,git_pos,i))
                 if (git < N_T) then
                    MSM_Indiv%a(1,git_pos+1,i) = rcx2(1)  !A prime
                    MSM_Indiv%a(3,git_pos+1,i) = rcx2(2)  !AIME prime
                 endif
              else ! not receiving ss yet
                 if ((git<70 .AND. pNRA<=70) .OR. (pNRA>N_T)) then  !starting from 70, everyone takes ss since no DRC in the SS policy.
                    !to calculate the value if not receiving SS
                    if (git < N_T) MSM_Indiv_Irecss%a(git_pos+1,i) = 0 !not receive ss at git
                    giRECSS_next = 1  !not take ss at git
                    call sub_ForwardInt(MSM_Indiv%a(4:6,git_pos,i),rcx2,MSM_Indiv%a(8,git_pos,i))
                    if (git < N_T) then
                       MSM_Indiv%a(1,git_pos+1,i) = rcx2(1)  !A prime
                       MSM_Indiv%a(3,git_pos+1,i) = rcx2(2)  !AIME prime
                    endif
                 endif
                 if (git >= pERA) then
                    !to calculate the value if receiving SS 
                    rcx3(3) = MSM_Indiv%a(6,git_pos,i) !L
                    giRECSS_next = 2  !take ss at git
                    call sub_ForwardInt(rcx3,rcx2,ev) !receive ss
                    if ( (ev>=MSM_Indiv%a(8,git_pos,i)) .OR. (git>=70 .AND. pNRA<=70) ) then
                       MSM_Indiv%a(4,git_pos,i) = rcx3(1) !c
                       MSM_Indiv%a(5,git_pos,i) = rcx3(2) !I
                       MSM_Indiv%a(8,git_pos,i) = ev
                       if (git < N_T) then
                          MSM_Indiv_Irecss%a(git_pos+1,i) = 1  !RECSS 
                          MSM_Indiv%a(1,git_pos+1,i) = rcx2(1)  !A prime
                          MSM_Indiv%a(3,git_pos+1,i) = rcx2(2)  !AIME prime
                       endif
                    endif
                 endif !if (git >= pERA)
              endif !if (giRECSS == 2)
          endif  !if (gcI_CF == -12)              

          !1-A,2-H,3-AIME,4-c,5-I,6-L,7-Labor,8-Value,9-Leisure error,10-XI error
          MSM_Indiv%a(7,git_pos,i) = max(0.0d0, 1.0d0-MSM_Indiv%a(5,git_pos,i)-MSM_Indiv%a(6,git_pos,i))  !Labor
          
          if (gcL_IsOutput) then
              if (abs(MSM_Indiv%a(6,git_pos,i)-1.0d0) < 0.1d0) then       !leisure==1
                  rtemp = exp(gr_leisure_taste(git,giHEALTH,giFS)+MSM_Indiv%a(9,git_pos,i))
              elseif (abs(MSM_Indiv%a(6,git_pos,i)-0.5d0) < 0.1d0) then   !leisure==0.5
                  rtemp = exp(gr_leisure_taste(git,giHEALTH,giFS)+MSM_Indiv%a(9,git_pos,i)) * plpt_t(git,giFS)
              else  !leisure==0
                  rtemp = 0.0d0
              endif
              MSM_Indiv%a(8,git_pos,i) = MSM_Indiv%a(8,git_pos,i) + rtemp !value including leisure
          endif   !if (gcL_IsOutput)
                    
          !state (giRECSS_next, Hprime, Aprime+spousal income next period) next period
          if ( git < N_T ) then
             
             giRECSS_next = MSM_Indiv_Irecss%a(git_pos+1,i) + 1
             
             if (MSM_Indiv%a(6,git_pos,i)>=0.95d0 .AND. (gcI_robust .NE. 8) .AND. (gcI_robust .NE. 10)  &
                        &  .AND. (gcI_robust .NE. 12)) then !not working no shock, except exo case
                rtemp = (1.0d0 - psigma_offjob)*MSM_Indiv%a(2,git_pos,i)  !not working
             else !working or exo case, there is shock
                if (gcL_XI_discrete) then
                   rtemp = 0.0d0 
                   do j = 1, N_XI_int
                      rtemp = rtemp + gr_XI_xw(j,2)
                      if (MSM_Indiv%a(10,git_pos,i) <= rtemp) then
                         rtemp2 = gr_XI_xw(j,3)
                         exit
                      endif
                   enddo ! do j=1, N_XI_int
                   if (j > N_XI_int) rtemp2 = gr_XI_xw(N_XI_int,3)
                else
                   rtemp2 = exp(MSM_Indiv%a(10,git_pos,i)*gr_XI_sd+gr_XI_mean) 
                endif   !if (gcL_XI_discrete)
                
                if (MSM_Indiv%a(6,git_pos,i)>=0.95d0) then  !not working
                    rtemp = psigma_offjob
                else
                    rtemp = psigma_onjob
                endif
                
                if (gcI_robust == 8 .OR. gcI_robust == 9) then  !exo (regardless of working) or lbd (working)
                   rtemp = (1.0d0 - rtemp)*MSM_Indiv%a(2,git_pos,i)+ &
                        & max(0.0d0, rtemp2*ppi*(1.0d0+palpha_I*dble(git-N_0)/gr_LBD_t+palpha_H*dble((git-N_0)**2)/gr_LBD_tsq))
                elseif (gcI_robust == 10 .OR. gcI_robust == 11) then  !exo v2 (regardless of working) lbd v2 (working)
                   rtemp = (1.0d0 - rtemp)*MSM_Indiv%a(2,git_pos,i)+ &
                        & max(0.0d0, rtemp2*ppi*(1.0d0+palpha_I*MSM_Indiv%a(2,git_pos,i)+palpha_H*MSM_Indiv%a(2,git_pos,i)**2 ))
                elseif (gcI_robust == 12 .OR. gcI_robust == 13) then  !exo Hone (regardless of working) lbd Hone (working)
                   rtemp = (1.0d0 - rtemp)*MSM_Indiv%a(2,git_pos,i)+ &
                        & max(0.0d0, rtemp2*ppi*(MSM_Indiv%a(2,git_pos,i)**palpha_H))   
                else  !working
                    rtemp = (1.0d0 - psigma_onjob)*MSM_Indiv%a(2,git_pos,i)+  &
                            & rtemp2*ppi*(MSM_Indiv%a(5,git_pos,i)**palpha_I)*(MSM_Indiv%a(2,git_pos,i)**palpha_H)  !pi * (I^alpha) * (H^beta)
                endif
             endif
             if (gcL_Iscapped) rtemp = min(sgridH%a(git_pos+1, N_H), rtemp)
             MSM_Indiv%a(2,git_pos+1,i) = rtemp  !Hprime
             !health status at git+1 is already set at SUBROUTINE INIT_HEALTH
             MSM_Indiv%a(1,git_pos+1,i) = MSM_Indiv%a(1,git_pos+1,i) + MSM_Indiv%a(16,git_pos+1,i)  !plus spouse income of next period
          else  !git>=N_T
             giRECSS_next = 2 
          endif  !if ( git < N_T ) then
          
          if (gcL_IsOutput) call sub_calbudget(i)
          
       enddo !do git = N_0, N_T
      
    enddo !i
    
    !==============================================================================================
    ! Calculate moments
    !==============================================================================================
    if (myid_mpi .NE. master_mpi) rtemp = SMM_MomentsDistance()
    
    return
    end subroutine SMM_SimulateProfiles

  !========================================================================================================
  ! Simulate lifecycle profiles all by interpolation, called by all nodes (master, head slaves and slaves)
  !gcL_CF_IW_decompose=.TRUE.
  !gcI_CF = 1...12
  !gcI_CF_IW=1....7
  !========================================================================================================
  subroutine SMM_SimulateProfiles_IW()
    implicit none
    integer(4) :: i, j, ix(3), ierr
    real(8) :: rtemp, rtemp2, rcx3(3), rcx2(2), ev, rVy(3), rptx, rVxx(2,2,2), rVxx2(2,2,2), rVxx3(2,2,2)
    
    !MSM_Indiv
    !1-A,2-H,3-AIME,4-c,5-I,6-L,7-Labor,8-Value,9-leisure error,10-XI error,
    !11-HEALTH shock,12-V0,13-V1, 14-Prob of leisure=0 or Phi(epsilon_tmin) in PT, 
    !15-FS error, 16-FS spouse inc
    !17-produced output, 18-ss tax, 19-medicare tax, 20-income tax, 21-received ss benefit, 22-ssdiinc
    !23-V(-1), 24-Phi(epsilon_tmax)
        
    !We need to start from N_0 where AIME=0
    MSM_Indiv%a(8,:,:) = - NR_BIG2 !VALUE
    MSM_Indiv%a(17:22,:,:) = 0.0d0  !17-produced output, 18-ss tax, 19-medicare tax, 20-income tax, 21-received ss benefit, 22-ssdiinc
    MSM_Indiv_Irecss%a(:,:) = 0    !Reset RECSS
    !do i = 1, N_Sim
    do i = gI_Indiv_Tasks_mpi%a(myid_mpi+1,1), gI_Indiv_Tasks_mpi%a(myid_mpi+1,2) 
       !========= Initial conditions ===============
       git = N_0
       git_pos = git - N_0 + 1
       MSM_Indiv%a(1,git_pos,i) = pinitA + MSM_Indiv_initX%a(i,1)      !initial Asset
       !MSM_Indiv%a(1,git_pos,i) = max(sgridA%a(git_pos,1), min(sgridA%a(git_pos,VN_A(git)), MSM_Indiv%a(1,git_pos,i))) !initial A
          
       MSM_Indiv%a(2,git_pos,i) = pinitHmean + pinitHsd * MSM_Indiv_initX%a(i,2) !initial H
       
       if (N_het_H0>1) then
          MSM_Indiv%a(2,git_pos,i) = MSM_Indiv%a(2,git_pos,i) + pinitH_a0*pa0 + pinitH_pi*ppi
       endif
       
       MSM_Indiv%a(2,git_pos,i) = exp(MSM_Indiv%a(2,git_pos,i))   !initial H
       
       MSM_Indiv%a(2,git_pos,i) = max(sgridH%a(git_pos,1)+1.0d-6, min(sgridH%a(git_pos,N_H)-1.0d-6, MSM_Indiv%a(2,git_pos,i))) !H
       MSM_Indiv%a(3,git_pos,i) = min(sgridAIME%a(git_pos,1)+1.0d-6, sgridAIME%a(git_pos,VN_AIME(git)))  !initial AIME at N_0
       
       !1 RECSS 2 HEALTH
       MSM_Indiv_Irecss%a(git_pos,i) = 0  !this is the status of current period: not receive ss at N_0
       
       !health status at N_0 is set at SUBROUTINE INIT_HEALTH
       !END of ========= Initial conditions ===============
       
       do git = N_0, N_T
          git_pos = git - N_0 + 1 
          call xPosadjust(git, MSM_Indiv%a(1:3,git_pos,i), ix)
          gA = MSM_Indiv%a(1,git_pos,i)
          gH = MSM_Indiv%a(2,git_pos,i)
          gAIME = MSM_Indiv%a(3,git_pos,i)
          giA = ix(1)
          giH = ix(2)
          giAIME = ix(3)
          giRECSS = MSM_Indiv_Irecss%a(git_pos,i) + 1   !index, 1-not RECSS, 2-RECSS
          if (giRECSS<1 .OR. giRECSS>2) then
             write(*,*) 'WRONG 1: i=',i,',git=',git,',giRECSS=',giRECSS, &
                      & ',MSM_Indiv_Irecss%a(N_0:N_T,i)=',MSM_Indiv_Irecss%a(N_0-N_0+1:N_T-N_0+1,i)
          endif
          giHEALTH = MSM_Indiv_Ihealth%a(git_pos,i)      !health status, 1-excellent, 2-good, 3-bad 4-disabled
          giFS = MSM_Indiv_IFS%a(git_pos,i)   !Family status: 1-3
          
          if (gcI_CF == -12 .AND. git == gcI_elas_t) then  !!unanticipated elasticity calculation. Use optimization when shock happens
              call SMM_SimulateProfiles_optimization(i)
          else   !use interpolation
              !First, get depsilon to decide leisure being 0 or 1
              if (Ipt == 0)  then
                 if (gcI_CF_IW==2 .OR. gcI_CF_IW==5 .OR. gcI_CF_IW==6 .OR. gcI_CF_IW==7) then
                     rVxx = depsilon0(0)%a(ix(1):(ix(1)+1),ix(3):(ix(3)+1),giHEALTH,giRECSS,giFS,ix(2):(ix(2)+1),git_pos)
                 else
                     rVxx = depsilon(0)%a(ix(1):(ix(1)+1),ix(3):(ix(3)+1),giHEALTH,giRECSS,giFS,ix(2):(ix(2)+1),git_pos)
                 endif
                 
                 if (gcL_triangle) then
                    rtemp = linear_inte_tetrahedra3(sgridA%a(git_pos,ix(1):(ix(1)+1)), &
                         &  sgridAIME%a(git_pos,ix(3):(ix(3)+1)), sgridH%a(git_pos,ix(2):(ix(2)+1)), &
                         &  rVxx,  &
                         &  (/gA,gAIME,gH/) ) !threshold
                 else
                    rtemp = linear_inte3(gcI_inte_IL_log, (/-1,-1,1/), sgridA%a(git_pos,ix(1):(ix(1)+1)), &
                         &  sgridAIME%a(git_pos,ix(3):(ix(3)+1)), sgridH%a(git_pos,ix(2):(ix(2)+1)), &
                         &  rVxx,  &
                         &  (/gA,gAIME,gH/) ) !threshold
                 endif
              else   !(Ipt == 1)
                  !always use gcL_triangle==.true.
                  
                  !!!! calculate V0, V1, and Vp, and then derive two thresholds.
                  !!assuming there is always stochastic shock
                  do j = 0, Ipt+1
                     if (gcI_CF_IW==2 .OR. gcI_CF_IW==5 .OR. gcI_CF_IW==6 .OR. gcI_CF_IW==7) then
                         rVxx = dvV0(j)%a(ix(1):(ix(1)+1),ix(3):(ix(3)+1),giHEALTH,giRECSS,giFS,ix(2):(ix(2)+1),git_pos)
                     else
                         rVxx = dvV(j)%a(ix(1):(ix(1)+1),ix(3):(ix(3)+1),giHEALTH,giRECSS,giFS,ix(2):(ix(2)+1),git_pos)
                     endif
                     
                     rVy(1+j) = linear_inte_tetrahedra3(sgridA%a(git_pos,ix(1):(ix(1)+1)), &
                         &  sgridAIME%a(git_pos,ix(3):(ix(3)+1)), sgridH%a(git_pos,ix(2):(ix(2)+1)), &
                         &  rVxx,  &
                         &  (/gA,gAIME,gH/) ) !v0
                  enddo !do j
                  rcx3(1) = rVy(1)  !leisure==0
                  rcx3(2) = rVy(2)  !leisure==1
                  rptx = rVy(3)     !leisure==2 (pt)
                  
                  rVy(3) = gr_leisure_taste(git,giHEALTH,giFS)
                  rcx3(3) = rcx3(1) - rcx3(2)  !leisure 0 - leisure 1
                  rVy(1) = min((rcx3(1)-rptx)/max(NR_EPS,plpt_t(git,giFS)), rcx3(3))   !gamma_tmin. Condition A8. Leisure 0 - p
                  rVy(2) = max((rptx-rcx3(2))/max(NR_EPS,(1.0d0-plpt_t(git,giFS))), rcx3(3))  !gamma_tmax. Leisure p - 1
                  
                   
                  rtemp2 = (log(max(NR_EPS, rVy(1)))-rVy(3)-grUVmultiplier_ln)/pa_epsilon !epsilon_tmin. /pa_epsilon is to get a standard normal with sd=1
                  rtemp = (log(max(NR_EPS, rVy(2)))-rVy(3)-grUVmultiplier_ln)/pa_epsilon !epsilon_tmax
              endif  !if (Ipt == 0)
                            
              if (Ipt == 0)  then  !full time only
                  if (MSM_Indiv%a(9,git_pos,i) <= rtemp) then !leisure = 0
                      MSM_Indiv%a(6,git_pos,i) = 0.0d0
                  else
                      MSM_Indiv%a(6,git_pos,i) = 1.0d0
                  endif
                  MSM_Indiv%a(14,git_pos,i) = CDF_Normal(rtemp, 0.0d0, 1.0d0) !prob of leisure=0
              else   !Ipt == 1, part time
                  if (MSM_Indiv%a(9,git_pos,i) <= rtemp2) then !leisure = 0
                      MSM_Indiv%a(6,git_pos,i) = 0.0d0
                  elseif (MSM_Indiv%a(9,git_pos,i) <= rtemp) then !leisure = 0.5
                      MSM_Indiv%a(6,git_pos,i) = 0.5d0
                  else
                      MSM_Indiv%a(6,git_pos,i) = 1.0d0
                  endif
                  MSM_Indiv%a(14,git_pos,i) = CDF_Normal(rtemp2, 0.0d0, 1.0d0)  !prob of leisure=0
                  MSM_Indiv%a(24,git_pos,i) = CDF_Normal(rtemp, 0.0d0, 1.0d0)   !prob of leisure=pt or 0
              endif  !if (Ipt == 0)
              
              if (gcI_CF_IW==2 .OR. gcI_CF_IW==5 .OR. gcI_CF_IW==6 .OR. gcI_CF_IW==7) then
                  rVxx = dvV0(0)%a(ix(1):(ix(1)+1),ix(3):(ix(3)+1),giHEALTH,giRECSS,giFS,ix(2):(ix(2)+1),git_pos)
                  rVxx2 = dvV0(1)%a(ix(1):(ix(1)+1),ix(3):(ix(3)+1),giHEALTH,giRECSS,giFS,ix(2):(ix(2)+1),git_pos)
                  if (Ipt==1) rVxx3 = dvV0(2)%a(ix(1):(ix(1)+1),ix(3):(ix(3)+1),giHEALTH,giRECSS,giFS,ix(2):(ix(2)+1),git_pos)
              else
                  rVxx = dvV(0)%a(ix(1):(ix(1)+1),ix(3):(ix(3)+1),giHEALTH,giRECSS,giFS,ix(2):(ix(2)+1),git_pos)
                  rVxx2 = dvV(1)%a(ix(1):(ix(1)+1),ix(3):(ix(3)+1),giHEALTH,giRECSS,giFS,ix(2):(ix(2)+1),git_pos)
                  if (Ipt==1) rVxx3 = dvV(2)%a(ix(1):(ix(1)+1),ix(3):(ix(3)+1),giHEALTH,giRECSS,giFS,ix(2):(ix(2)+1),git_pos)
              endif
              
              if (gcL_triangle) then
                 MSM_Indiv%a(12,git_pos,i) = linear_inte_tetrahedra3(sgridA%a(git_pos,ix(1):(ix(1)+1)), &
                         &  sgridAIME%a(git_pos,ix(3):(ix(3)+1)), sgridH%a(git_pos,ix(2):(ix(2)+1)), &
                         &  rVxx,  &
                         &  (/gA,gAIME,gH/) ) !V0
                 MSM_Indiv%a(13,git_pos,i) = linear_inte_tetrahedra3(sgridA%a(git_pos,ix(1):(ix(1)+1)), &
                         &  sgridAIME%a(git_pos,ix(3):(ix(3)+1)), sgridH%a(git_pos,ix(2):(ix(2)+1)), &
                         &  rVxx2,  &
                         &  (/gA,gAIME,gH/) ) !V1
                 if (Ipt==1) then
                     MSM_Indiv%a(23,git_pos,i) = linear_inte_tetrahedra3(sgridA%a(git_pos,ix(1):(ix(1)+1)), &
                         &  sgridAIME%a(git_pos,ix(3):(ix(3)+1)), sgridH%a(git_pos,ix(2):(ix(2)+1)), &
                         &  rVxx3,  &
                         &  (/gA,gAIME,gH/) ) !V(2)
                 endif
              else
                 MSM_Indiv%a(12,git_pos,i) = linear_inte3(gcI_inte_v_log, (/1,1,1/), sgridA%a(git_pos,ix(1):(ix(1)+1)), &
                         &  sgridAIME%a(git_pos,ix(3):(ix(3)+1)), sgridH%a(git_pos,ix(2):(ix(2)+1)), &
                         &  rVxx,  &
                         &  (/gA,gAIME,gH/) ) !V0
                 MSM_Indiv%a(13,git_pos,i) = linear_inte3(gcI_inte_v_log, (/1,1,1/), sgridA%a(git_pos,ix(1):(ix(1)+1)), &
                         &  sgridAIME%a(git_pos,ix(3):(ix(3)+1)), sgridH%a(git_pos,ix(2):(ix(2)+1)), &
                         &  rVxx2,  &
                         &  (/gA,gAIME,gH/) ) !V1
                 if (Ipt==1) then
                     MSM_Indiv%a(23,git_pos,i) = linear_inte3(gcI_inte_v_log, (/1,1,1/), sgridA%a(git_pos,ix(1):(ix(1)+1)), &
                         &  sgridAIME%a(git_pos,ix(3):(ix(3)+1)), sgridH%a(git_pos,ix(2):(ix(2)+1)), &
                         &  rVxx3,  &
                         &  (/gA,gAIME,gH/) ) !V(2)
                 endif
              endif 
              
              if (giRECSS == 2) then  !already receiving ss 
                 !the filter makes sure that now git >= pERA
                 if (git < N_T) MSM_Indiv_Irecss%a(git_pos+1,i) = 1  !receiving ss is absorbing
                 giRECSS_next = 2
                 call sub_ForwardInt_IW(MSM_Indiv%a(4:6,git_pos,i),rcx2,MSM_Indiv%a(8,git_pos,i))
                 if (git < N_T) then
                    MSM_Indiv%a(1,git_pos+1,i) = rcx2(1)  !A prime
                    MSM_Indiv%a(3,git_pos+1,i) = rcx2(2)  !AIME prime
                 endif
              else ! not receiving ss yet
                 if ((git<70 .AND. pNRA<=70) .OR. (pNRA>N_T)) then  !starting from 70, everyone takes ss since no DRC in the SS policy.
                    !to calculate the value if not receiving SS
                    if (git < N_T) MSM_Indiv_Irecss%a(git_pos+1,i) = 0 !not receive ss at git
                    giRECSS_next = 1  !not take ss at git
                    call sub_ForwardInt_IW(MSM_Indiv%a(4:6,git_pos,i),rcx2,MSM_Indiv%a(8,git_pos,i))
                    if (git < N_T) then
                       MSM_Indiv%a(1,git_pos+1,i) = rcx2(1)  !A prime
                       MSM_Indiv%a(3,git_pos+1,i) = rcx2(2)  !AIME prime
                    endif
                 endif
                 if (git >= pERA) then
                    !to calculate the value if receiving SS 
                    rcx3(3) = MSM_Indiv%a(6,git_pos,i) !L
                    giRECSS_next = 2  !take ss at git
                    call sub_ForwardInt_IW(rcx3,rcx2,ev) !receive ss
                    if ( (ev>=MSM_Indiv%a(8,git_pos,i)) .OR. (git>=70 .AND. pNRA<=70) ) then
                       MSM_Indiv%a(4,git_pos,i) = rcx3(1) !c
                       MSM_Indiv%a(5,git_pos,i) = rcx3(2) !I
                       MSM_Indiv%a(8,git_pos,i) = ev
                       if (git < N_T) then
                          MSM_Indiv_Irecss%a(git_pos+1,i) = 1  !RECSS 
                          MSM_Indiv%a(1,git_pos+1,i) = rcx2(1)  !A prime
                          MSM_Indiv%a(3,git_pos+1,i) = rcx2(2)  !AIME prime
                       endif
                    endif
                 endif !if (git >= pERA)
              endif !if (giRECSS == 2)
          endif  !if (gcI_CF == -12)              

          !1-A,2-H,3-AIME,4-c,5-I,6-L,7-Labor,8-Value,9-Leisure error,10-XI error
          MSM_Indiv%a(7,git_pos,i) = max(0.0d0, 1.0d0-MSM_Indiv%a(5,git_pos,i)-MSM_Indiv%a(6,git_pos,i))  !Labor
          
          if (gcL_IsOutput) then
              if (abs(MSM_Indiv%a(6,git_pos,i)-1.0d0) < 0.1d0) then       !leisure==1
                  rtemp = exp(gr_leisure_taste(git,giHEALTH,giFS)+MSM_Indiv%a(9,git_pos,i))
              elseif (abs(MSM_Indiv%a(6,git_pos,i)-0.5d0) < 0.1d0) then   !leisure==0.5
                  rtemp = exp(gr_leisure_taste(git,giHEALTH,giFS)+MSM_Indiv%a(9,git_pos,i)) * plpt_t(git,giFS)
              else  !leisure==0
                  rtemp = 0.0d0
              endif
              MSM_Indiv%a(8,git_pos,i) = MSM_Indiv%a(8,git_pos,i) + rtemp !value including leisure
          endif   !if (gcL_IsOutput)
          
          
          !state (giRECSS_next, Hprime, Aprime+spousal income next period) next period
          if ( git < N_T ) then
             
             giRECSS_next = MSM_Indiv_Irecss%a(git_pos+1,i) + 1
             
             if (MSM_Indiv%a(6,git_pos,i)>=0.95d0 .AND. (gcI_robust .NE. 8) .AND. (gcI_robust .NE. 10)  &
                        &  .AND. (gcI_robust .NE. 12)) then !not working no shock, except exo case
                rtemp = (1.0d0 - psigma_offjob)*MSM_Indiv%a(2,git_pos,i)  !not working
             else !working or exo case, there is shock
                if (gcL_XI_discrete) then
                   rtemp = 0.0d0 
                   do j = 1, N_XI_int
                      rtemp = rtemp + gr_XI_xw(j,2)
                      if (MSM_Indiv%a(10,git_pos,i) <= rtemp) then
                         rtemp2 = gr_XI_xw(j,3)
                         exit
                      endif
                   enddo ! do j=1, N_XI_int
                   if (j > N_XI_int) rtemp2 = gr_XI_xw(N_XI_int,3)
                else
                   rtemp2 = exp(MSM_Indiv%a(10,git_pos,i)*gr_XI_sd+gr_XI_mean) 
                endif   !if (gcL_XI_discrete)
                
                if (MSM_Indiv%a(6,git_pos,i)>=0.95d0) then  !not working
                    rtemp = psigma_offjob
                else
                    rtemp = psigma_onjob
                endif
                
                if (gcI_robust == 8 .OR. gcI_robust == 9) then  !exo (regardless of working) or lbd (working)
                   rtemp = (1.0d0 - rtemp)*MSM_Indiv%a(2,git_pos,i)+ &
                        & max(0.0d0, rtemp2*ppi*(1.0d0+palpha_I*dble(git-N_0)/gr_LBD_t+palpha_H*dble((git-N_0)**2)/gr_LBD_tsq))
                elseif (gcI_robust == 10 .OR. gcI_robust == 11) then  !exo v2 (regardless of working) lbd v2 (working)
                   rtemp = (1.0d0 - rtemp)*MSM_Indiv%a(2,git_pos,i)+ &
                        & max(0.0d0, rtemp2*ppi*(1.0d0+palpha_I*MSM_Indiv%a(2,git_pos,i)+palpha_H*MSM_Indiv%a(2,git_pos,i)**2 ))
                elseif (gcI_robust == 12 .OR. gcI_robust == 13) then  !exo Hone (regardless of working) lbd Hone (working)
                   rtemp = (1.0d0 - rtemp)*MSM_Indiv%a(2,git_pos,i)+ &
                        & max(0.0d0, rtemp2*ppi*(MSM_Indiv%a(2,git_pos,i)**palpha_H))   
                else  !working
                    rtemp = (1.0d0 - psigma_onjob)*MSM_Indiv%a(2,git_pos,i)+  &
                            & rtemp2*ppi*(MSM_Indiv%a(5,git_pos,i)**palpha_I)*(MSM_Indiv%a(2,git_pos,i)**palpha_H)  !pi * (I^alpha) * (H^beta)
                endif
             endif
             if (gcL_Iscapped) rtemp = min(sgridH%a(git_pos+1, N_H), rtemp)
             MSM_Indiv%a(2,git_pos+1,i) = rtemp  !Hprime
             !health status at git+1 is already set at SUBROUTINE INIT_HEALTH
             MSM_Indiv%a(1,git_pos+1,i) = MSM_Indiv%a(1,git_pos+1,i) + MSM_Indiv%a(16,git_pos+1,i)  !plus spouse income of next period
          else  !git>=N_T
             giRECSS_next = 2 
          endif  !if ( git < N_T ) then
          
          if (gcL_IsOutput) call sub_calbudget(i)
          
       enddo !do git = N_0, N_T
      
    enddo !i
    
    !==============================================================================================
    ! Calculate moments
    !==============================================================================================
    if (myid_mpi .NE. master_mpi) rtemp = SMM_MomentsDistance()
    
    return
 end subroutine SMM_SimulateProfiles_IW

     
  !========================================================================================================
  ! Simulate lifecycle profiles all by optimization, called by all nodes (master, head slaves and slaves)
  !========================================================================================================
  subroutine SMM_SimulateProfiles_optimization(i)
    implicit none
    INTEGER(4), INTENT(IN) :: i

    integer(4) :: iRecss(0:1), ierr
    real(8) :: rtemp, rcx2(2), ev, rc(0:1), rI, rcx(3), rVy(2), repsilon
    
    !MSM_Indiv
    !1-A,2-H,3-AIME,4-c,5-I,6-L,7-Labor,8-Value,9-leisure error,10-XI error
    !11-HEALTH shock,12-V0,13-V1, 14-Prob of leisure=0, 
    !15-FS error, 16-FS spouse inc
    !17-produced output, 18-payroll tax, 19-income tax, 20-received ss benefit
    
    !calculate the values of 0 and 1 by optimization
    
    !======================the innerest loop====================================
    !choose (c, I, L, recss) to maximize value function
    !===========================================================================
    if (giRECSS == 2) then  !already receiving ss.
       giRECSS_next = 2  !receiving ss is absorbing. the filter makes sure that now git >= pERA
       if (git < N_T) MSM_Indiv_Irecss%a(git_pos+1,i) = 1  !receiving ss is absorbing
       call dcmaxee_getVAL(MSM_Indiv%a(12,git_pos,i), rc(0), rI, 0)    !case with leisure = 0
       call dcmaxee_getVAL(MSM_Indiv%a(13,git_pos,i), rc(1), rtemp, 1)  !case with leisure = 1
    else ! not receiving ss yet
       if (git<70) then
          !===== case with leisure = 0
          giRECSS_next = 1 !not receive ss at git
          call dcmaxee_getVAL(MSM_Indiv%a(12,git_pos,i), rc(0), rI, 0)    !case with leisure = 0
          iRecss(0) = giRECSS_next-1
          if (git >= pERA) then ! check if they should start collecting social security benefit at this period
             giRECSS_next = 2
             ev = - NR_BIG2
             call dcmaxee_getVAL(ev, rcx2(1), rcx2(2), 0)
             if (ev >= MSM_Indiv%a(12,git_pos,i)) then
                rc(0) = rcx2(1)
                rI = rcx2(2)
                iRecss(0) = giRECSS_next-1
                MSM_Indiv%a(12,git_pos,i) = ev
             endif
          endif !if (git >= pERA)
          
          !=====case with leisure = 1
          giRECSS_next = 1 !not receive ss at git
          call dcmaxee_getVAL(MSM_Indiv%a(13,git_pos,i), rc(1), rtemp, 1)
          iRecss(1) = giRECSS_next-1
          if (git >= pERA) then
             giRECSS_next = 2
             ev = - NR_BIG2
             call dcmaxee_getVAL(ev, rcx2(1), rcx2(2), 1)
             if (ev >= MSM_Indiv%a(13,git_pos,i)) then
                rc(1) = rcx2(1)
                iRecss(1) = giRECSS_next-1
                MSM_Indiv%a(13,git_pos,i) = ev
             endif
          endif !if (git >= pERA)
       else  !git>=70. Always take ss
          !===== case with leisure = 0
          giRECSS_next = 2 !always take ss
          call dcmaxee_getVAL(MSM_Indiv%a(12,git_pos,i), rc(0), rI, 0)
          iRecss(0) = giRECSS_next-1
          !=====case with leisure = 1
          giRECSS_next = 2 !always take ss
          call dcmaxee_getVAL(MSM_Indiv%a(13,git_pos,i), rc(1), rtemp, 1)
          iRecss(1) = giRECSS_next-1
       endif !if (git<70)
    endif !if (giRECSS == 2)
    !========================End of the innerest loop=================================
    !calculate Emax and epsilon
    rVy(2) = gr_leisure_taste(git,giHEALTH,giFS)
                          
    rcx(1) = MSM_Indiv%a(12,git_pos,i)   !leisure=0
    rcx(2) = MSM_Indiv%a(13,git_pos,i)   !leisure=1
    if (abs(pa_epsilon) <= 0.0d0) then  !no stochastic shock
        rcx(3) = exp(rVy(2))
        rcx(2) = rcx(2) + rcx(3)
        if (rcx(2) >= rcx(1)) then
            rVy(1) = 0.0d0  !prob of (leisure=0)
            repsilon = - 1.0d3 !leisure=1 for sure
        else
            rVy(1) = 1.0d0  !prob of (leisure=0)
            repsilon = 1.0d3  !leisure=0 for sure
        endif
        ev = rVy(1) * rcx(1) + (1.0d0-rVy(1)) * rcx(2) 
    else   !if (pa_epsilon>0.0d0) then
        rcx(3) = rcx(1) - rcx(2)  !leisure 0 - leisure 1
        if (rcx(3) < NR_EPS) rcx(3) = NR_EPS
        rcx(3) = (log(rcx(3))-rVy(2)-grUVmultiplier_ln)/pa_epsilon  !depsilon
        repsilon = rcx(3)
        rVy(1) = CDF_Normal(rcx(3), 0.0d0, 1.0d0) !Phi(depsilon)
        rcx(3) = gr_leisure_condexp(git,giHEALTH,giFS) * CDF_Normal(pa_epsilon - rcx(3), 0.0d0, 1.0d0)                    
        ev = rVy(1) * rcx(1) + (1.0d0 - rVy(1)) * rcx(2) + rcx(3)
    endif  !if (abs(pa_epsilon)<NR_EPS) then
    MSM_Indiv%a(14,git_pos,i) = CDF_Normal(repsilon, 0.0d0, 1.0d0)   !14-Prob of leisure=0
    
    !Now, decide leisure being 0 or 1  
    if (MSM_Indiv%a(9,git_pos,i) <= repsilon) then !leisure = 0
        MSM_Indiv%a(6,git_pos,i) = 0.0d0
        MSM_Indiv%a(4,git_pos,i) = rc(0)
        MSM_Indiv%a(5,git_pos,i) = rI
        giRECSS_next = iRecss(0) + 1
    else
        MSM_Indiv%a(6,git_pos,i) = 1.0d0
        MSM_Indiv%a(4,git_pos,i) = rc(1)
        MSM_Indiv%a(5,git_pos,i) = 0.0d0
        giRECSS_next = iRecss(1) + 1
    endif
    !1-A,2-H,3-AIME,4-c,5-I,6-L,7-Labor,8-Value,9-leisure error,10-XI error,
    !11-HEALTH shock,12-V0,13-V1, 14-prob of leisure=0
    !The next period state variables
    rtemp = pw(git) * gH * (1.0d0 - MSM_Indiv%a(5,git_pos,i)) * (1.0d0 - MSM_Indiv%a(6,git_pos,i))      !Labor income
    call A_AIME_prime(0.0d0, rtemp, rcx2)  !given git, gA, gAIME, giRECSS, giRECSS_next (current policy) 
    rcx2(1) = rcx2(1) - MSM_Indiv%a(4,git_pos,i) !Aprime
    !It is gurranteed Aprime is no lower than the lower bound of asset
    if (rcx2(1) < sgridAnrlb%a(git_pos+1)) then
        MSM_Indiv%a(4,git_pos,i) = grCF   !consumption
        rcx2(1) = sgridAnrlb%a(git_pos+1)
    endif
    if (gcL_Iscapped .AND. (git < N_T)) then
        rcx2(1) = min(sgridA%a(git_pos+1, VN_A(git+1)), rcx2(1))
        rcx2(2) = min(sgridAIME%a(git_pos+1, VN_AIME(git+1)), rcx2(2)) !AIME prime
    endif
    if (git < N_T) then
        MSM_Indiv%a(1,git_pos+1,i) = rcx2(1)  !A prime
        MSM_Indiv%a(3,git_pos+1,i) = rcx2(2)  !AIME prime
    endif
    
    !Value
    MSM_Indiv%a(8,git_pos,i) = ufunc_c(MSM_Indiv%a(4,git_pos,i)) + cmaxee_ev(MSM_Indiv%a(5:6,git_pos,i), rcx2)
    
    if (gcL_IsOutput) then
        if (abs(MSM_Indiv%a(6,git_pos,i)-1.0d0) < 0.1d0) then       !leisure==1
            rtemp = exp(gr_leisure_taste(git,giHEALTH,giFS)+MSM_Indiv%a(9,git_pos,i))
        elseif (abs(MSM_Indiv%a(6,git_pos,i)-0.5d0) < 0.1d0) then   !leisure==0.5
            rtemp = exp(gr_leisure_taste(git,giHEALTH,giFS)+MSM_Indiv%a(9,git_pos,i)) * plpt_t(git,giFS)
        else  !leisure==0
            rtemp = 0.0d0
        endif
        MSM_Indiv%a(8,git_pos,i) = MSM_Indiv%a(8,git_pos,i) + rtemp !value including leisure
    endif   !if (gcL_IsOutput)
    
    return
    end subroutine SMM_SimulateProfiles_optimization

  !==========================================================================================
  !when calculating the forward life-cycle path, do the interpolation within the period
  ! be careful here that "irecss_" is the state variable while "policy_ret_" is the policy
  ! The rnext_(1) is the asset at the end of Period git
  ! rnext_(2) is the AIME next period
  ! rv_ includes the utility from consumption, SSA and continuation (NO Leisure part)
  !==========================================================================================
  subroutine sub_ForwardInt(rc_, rnext_, rv_)
     implicit none
     real(8), dimension(:), intent(inout) :: rc_     !C,I,L
     real(8), dimension(:), intent(out) :: rnext_    !(A,AIME) prime
     real(8), intent(out) :: rv_
     integer(4) :: ix(3), i, j, iLeisure
     real(8) :: sx(3), rtemp

     if (rc_(3)<0.3d0) then
        iLeisure = 0
     elseif (rc_(3)>0.8d0) then
        iLeisure = 1 
     else
        iLeisure = 2 
     endif
     
     sx = (/gA,gH,gAIME/)   !state variables
     call xPosadjust(git, sx, ix)
     
     !Investment
     if (iLeisure==1 .OR. (gcI_robust>=8 .AND. gcI_robust<=13)) then  !L
        rc_(2) = 0.0d0    !I
     else
        if (gcL_triangle) then 
           rc_(2) = linear_inte_tetrahedra3(sgridA%a(git_pos,ix(1):(ix(1)+1)), &
                           &  sgridAIME%a(git_pos,ix(3):(ix(3)+1)), sgridH%a(git_pos,ix(2):(ix(2)+1)),  &
                           &  dcI(iLeisure/2)%a(ix(1):(ix(1)+1),ix(3):(ix(3)+1),giHEALTH,giRECSS,giFS, &
                           &  ix(2):(ix(2)+1),git_pos),(/sx(1),sx(3),sx(2)/)) !I 
        else
           rc_(2) = linear_inte3(gcI_inte_IL_log, (/0,0,0/), sgridA%a(git_pos,ix(1):(ix(1)+1)), &
                           &  sgridAIME%a(git_pos,ix(3):(ix(3)+1)), sgridH%a(git_pos,ix(2):(ix(2)+1)),  &
                           &  dcI(iLeisure/2)%a(ix(1):(ix(1)+1),ix(3):(ix(3)+1),giHEALTH,giRECSS,giFS,  &
                           &  ix(2):(ix(2)+1),git_pos),(/sx(1),sx(3),sx(2)/)) !I
        endif
        if ((rc_(2)>gr_p_Ibound(iLeisure,2)) .or. (rc_(2)<0.0d0)) then
           rc_(2) = max(0.0d0, min(gr_p_Ibound(iLeisure,2), rc_(2)))  !I
        endif
     endif
     
     !Consumption
     if (gcL_triangle) then
        rc_(1) = linear_inte_tetrahedra3(sgridA%a(git_pos,ix(1):(ix(1)+1)), &
            &  sgridAIME%a(git_pos,ix(3):(ix(3)+1)), sgridH%a(git_pos,ix(2):(ix(2)+1)), &
            &  dcc(iLeisure)%a(ix(1):(ix(1)+1),ix(3):(ix(3)+1),giHEALTH,giRECSS,giFS, &
            &  ix(2):(ix(2)+1),git_pos),(/sx(1),sx(3),sx(2)/)) !c
                           !&  dcc(iLeisure)%a(git_pos,ix(1):(ix(1)+1),ix(2):(ix(2)+1),ix(3):(ix(3)+1),irecss_+1),sx) !c
     else
        rc_(1) = linear_inte3(gcI_inte_c_log, (/1,1,1/), sgridA%a(git_pos,ix(1):(ix(1)+1)), &
            &  sgridAIME%a(git_pos,ix(3):(ix(3)+1)), sgridH%a(git_pos,ix(2):(ix(2)+1)), &
            &  dcc(iLeisure)%a(ix(1):(ix(1)+1),ix(3):(ix(3)+1),giHEALTH,giRECSS,giFS, &
            &  ix(2):(ix(2)+1),git_pos),(/sx(1),sx(3),sx(2)/)) !c
     endif
     
     if (rc_(1) <= 0.0d0) then
         rc_(1) = linear_inte3(1, (/0,0,0/), sgridA%a(git_pos,ix(1):(ix(1)+1)), &
               &  sgridAIME%a(git_pos,ix(3):(ix(3)+1)), sgridH%a(git_pos,ix(2):(ix(2)+1)), &
               &  dcc(iLeisure)%a(ix(1):(ix(1)+1),ix(3):(ix(3)+1),giHEALTH,giRECSS,giFS, &
               &  ix(2):(ix(2)+1),git_pos),(/sx(1),sx(3),sx(2)/)) !c
                         !&  dcc(iLeisure)%a(git_pos,ix(1):(ix(1)+1),ix(2):(ix(2)+1),ix(3):(ix(3)+1),irecss_+1),sx) !c
     endif

     !The next period state variables
     if (iLeisure==1) then  !Leisure=1
        rtemp = 0.0d0      !Labor income
     elseif (iLeisure==0) then  !Leisure=0, full time
        rtemp = pw(git) * sx(2) * (1.0d0 - rc_(2))      !Labor income 
     else  !leisure=2, part time
        rtemp = pw(git) * sx(2) * (0.5d0 - rc_(2))      !Labor income  
     endif
     
     call A_AIME_prime(0.0d0, rtemp, rnext_)  !given git, gA, gAIME, giRECSS, giRECSS_next (current policy) 
     
     rnext_(1) = rnext_(1) - rc_(1) !Aprime     
     !It is gurranteed Aprime is no lower than the lower bound of asset
     if (rnext_(1) < sgridAnrlb%a(git_pos+1)) then
        rc_(1) = grCF   !consumption
        rnext_(1) = sgridAnrlb%a(git_pos+1)
     endif
     
     if (gcL_Iscapped .AND. (git < N_T)) then
        rnext_(1) = min(sgridA%a(git_pos+1, VN_A(git+1)), rnext_(1))
        rnext_(2) = min(sgridAIME%a(git_pos+1, VN_AIME(git+1)), rnext_(2)) !AIME prime
     endif
    
     !Value
     rv_ = ufunc_c(rc_(1)) + cmaxee_ev(rc_(2:3), rnext_)
     
  end subroutine sub_ForwardInt
  
  !gcL_CF_IW_decompose=.TRUE.
  !gcI_CF = 1...27  
  !gcI_CF_IW=1...7
  subroutine sub_ForwardInt_IW(rc_, rnext_, rv_)
     implicit none
     real(8), dimension(:), intent(inout) :: rc_     !C,I,L
     real(8), dimension(:), intent(out) :: rnext_    !(A,AIME) prime
     real(8), intent(out) :: rv_
     integer(4) :: ix(3), i, j, iLeisure
     real(8) :: sx(3), rtemp, rVxx(2,2,2)

     if (rc_(3)<0.3d0) then
        iLeisure = 0
     elseif (rc_(3)>0.8d0) then
        iLeisure = 1 
     else
        iLeisure = 2 
     endif
     
     sx = (/gA,gH,gAIME/)   !state variables
     call xPosadjust(git, sx, ix)
     
     !Investment
     if (iLeisure==1 .OR. (gcI_robust>=8 .AND. gcI_robust<=13)) then  !L
        rc_(2) = 0.0d0    !I
     else
        if (gcI_CF_IW==1 .OR. gcI_CF_IW==4 .OR. gcI_CF_IW==6 .OR. gcI_CF_IW==7) then
            rVxx = dcI0(iLeisure/2)%a(ix(1):(ix(1)+1),ix(3):(ix(3)+1),giHEALTH,giRECSS,giFS,ix(2):(ix(2)+1),git_pos)
        else
            rVxx = dcI(iLeisure/2)%a(ix(1):(ix(1)+1),ix(3):(ix(3)+1),giHEALTH,giRECSS,giFS,ix(2):(ix(2)+1),git_pos)
        endif
        
        if (gcL_triangle) then 
           rc_(2) = linear_inte_tetrahedra3(sgridA%a(git_pos,ix(1):(ix(1)+1)), &
                           &  sgridAIME%a(git_pos,ix(3):(ix(3)+1)), sgridH%a(git_pos,ix(2):(ix(2)+1)),  &
                           &  rVxx,(/sx(1),sx(3),sx(2)/)) !I 
        else
           rc_(2) = linear_inte3(gcI_inte_IL_log, (/0,0,0/), sgridA%a(git_pos,ix(1):(ix(1)+1)), &
                           &  sgridAIME%a(git_pos,ix(3):(ix(3)+1)), sgridH%a(git_pos,ix(2):(ix(2)+1)),  &
                           &  rVxx,(/sx(1),sx(3),sx(2)/)) !I
        endif
        if ((rc_(2)>gr_p_Ibound(iLeisure,2)) .or. (rc_(2)<0.0d0)) then
           rc_(2) = max(0.0d0, min(gr_p_Ibound(iLeisure,2), rc_(2)))  !I
        endif
     endif
     
     !Consumption
     if (gcI_CF_IW==3 .OR. gcI_CF_IW==4 .OR. gcI_CF_IW==5 .OR. gcI_CF_IW==7) then
         rVxx = dcc0(iLeisure)%a(ix(1):(ix(1)+1),ix(3):(ix(3)+1),giHEALTH,giRECSS,giFS,ix(2):(ix(2)+1),git_pos)
     else
         rVxx = dcc(iLeisure)%a(ix(1):(ix(1)+1),ix(3):(ix(3)+1),giHEALTH,giRECSS,giFS,ix(2):(ix(2)+1),git_pos)
     endif
     
     if (gcL_triangle) then
        rc_(1) = linear_inte_tetrahedra3(sgridA%a(git_pos,ix(1):(ix(1)+1)), &
            &  sgridAIME%a(git_pos,ix(3):(ix(3)+1)), sgridH%a(git_pos,ix(2):(ix(2)+1)), &
            &  rVxx,(/sx(1),sx(3),sx(2)/)) !c
                           !&  dcc(iLeisure)%a(git_pos,ix(1):(ix(1)+1),ix(2):(ix(2)+1),ix(3):(ix(3)+1),irecss_+1),sx) !c
     else
        rc_(1) = linear_inte3(gcI_inte_c_log, (/1,1,1/), sgridA%a(git_pos,ix(1):(ix(1)+1)), &
            &  sgridAIME%a(git_pos,ix(3):(ix(3)+1)), sgridH%a(git_pos,ix(2):(ix(2)+1)), &
            &  rVxx,(/sx(1),sx(3),sx(2)/)) !c
     endif
     
     if (rc_(1) <= 0.0d0) then
         rc_(1) = linear_inte3(1, (/0,0,0/), sgridA%a(git_pos,ix(1):(ix(1)+1)), &
               &  sgridAIME%a(git_pos,ix(3):(ix(3)+1)), sgridH%a(git_pos,ix(2):(ix(2)+1)), &
               &  rVxx,(/sx(1),sx(3),sx(2)/)) !c
     endif

     !The next period state variables
     if (iLeisure==1) then  !Leisure=1
        rtemp = 0.0d0      !Labor income
     elseif (iLeisure==0) then  !Leisure=0, full time
        rtemp = pw(git) * sx(2) * (1.0d0 - rc_(2))      !Labor income 
     else  !leisure=2, part time
        rtemp = pw(git) * sx(2) * (0.5d0 - rc_(2))      !Labor income  
     endif
     
     call A_AIME_prime(0.0d0, rtemp, rnext_)  !given git, gA, gAIME, giRECSS, giRECSS_next (current policy) 
     
     rnext_(1) = rnext_(1) - rc_(1) !Aprime     
     !It is gurranteed Aprime is no lower than the lower bound of asset
     if (rnext_(1) < sgridAnrlb%a(git_pos+1)) then
        rc_(1) = grCF   !consumption
        rnext_(1) = sgridAnrlb%a(git_pos+1)
     endif
     
     if (gcL_Iscapped .AND. (git < N_T)) then
        rnext_(1) = min(sgridA%a(git_pos+1, VN_A(git+1)), rnext_(1))
        rnext_(2) = min(sgridAIME%a(git_pos+1, VN_AIME(git+1)), rnext_(2)) !AIME prime
     endif
    
     !Value
     rv_ = ufunc_c(rc_(1)) + cmaxee_ev(rc_(2:3), rnext_)
     
  end subroutine sub_ForwardInt_IW
  
  !this subroutine is called only when (gcL_IsOutput), in all nodes
  !It is to calculate !17-produced output=wage, 18-payroll tax, 19-income tax, 20-received ss benefit in MSM_Indiv for git
  !MSM_Indiv
  !1-A,2-H,3-AIME,4-c,5-I,6-L,7-Labor,8-Value,9-leisure error,10-XI error
  !11-HEALTH shock,12-V0,13-V1, 14-Prob of leisure=0, 
  !15-FS error, 16-FS spouse inc
  !17-produced output, 18-ss tax, 19-medicare tax, 20-income tax, 21-received ss benefit, 22-ssdiinc
  
  !similar as A_AIME_prime@cmaxee, but no need to update prime.
  subroutine sub_calbudget(i)
      implicit none
      integer(4), INTENT(IN) :: i
      integer(4) :: diffage
      real(8) :: rtax_ss, rtax_medicare, rtax_ss_cap, winc, ssinc, ssdiinc, pia, rothinc,  &
              &  rtinc, rtemp, rtemp2, rtemp3 
      
      if ( (gcI_CF==3) .OR. (gcI_CF==7) .OR. (gcI_CF==10) .OR. (gcI_CF==73)) then
          rtax_ss = 0.0d0
          rtax_medicare = 0.0d0
      else
          rtax_ss = 0.062d0
          rtax_medicare = 0.0145d0 
      endif
      rtax_ss_cap = 87900.0d0/grTotHours
    
      if (MSM_Indiv%a(6,git_pos,i)<0.9d0) then   !working
          winc = pw(git) * MSM_Indiv%a(2,git_pos,i)*MSM_Indiv%a(7,git_pos,i)         !actual labor income
      else
          winc = 0.0d0
      endif
      MSM_Indiv%a(17,git_pos,i) = winc  !17-produced output
      
      ssinc = 0.0d0
      ssdiinc = 0.0d0
      pia = AIME2PIA(gAIME)  !Effective PIA and  effective AIME   
      rothinc = 0.0d0
      
      !====================================================================
      ! Social security rules:
      ! 1. early retirement
      ! 2. delayed retirement credit 
      ! output: ssinc, 
      !         current_AIMEprime_adj (use this to adjust AIME next period)
      !=====================================================================
      if (giHEALTH==4 .AND. git<pNRA .AND. winc<gr_SSDI_SGA) then  !SSDI and SSI
          ssdiinc = pia * 12.0d0 !translate monthly PIA to yearly SSDI income
          ssdiinc = max(ssdiinc, gr_SSI)
      elseif (giRECSS_next == 2) then  !receiving ss
          if (giRECSS == 2) then
             ! has already taken ss in previous period
             ssinc = pia * 12.0d0 !translate monthly PIA to yearly ss income
          elseif (giRECSS == 1) then
             ! didn't take ss in previous period
             ! current period is the first time to take ss
             if (git < pERA) then
                 ssinc = 0.0d0
             elseif (git < pNRA) then  !early retirement
                 diffage = pNRA - git
                 if (diffage <= 3) then
                     rtemp = pia * (1.0d0 - 0.0667d0 * dble(diffage))
                 else
                     rtemp = pia * (0.8d0 - 0.05d0 * dble(diffage - 3))
                 endif
                 ssinc = rtemp * 12.0d0 !translate monthly PIA to yearly ss income
             elseif ( gcL_IsDRC .AND. (git > pNRA) ) then  !DRC: delayed retirement credit
                 rtemp = 1.0d0+dble(max(0, min(git,69)-pNRA))*0.06d0
                 ssinc = pia * rtemp * 12.0d0 !translate monthly PIA to yearly ss income
             else  !retire at NRA or (after NRA and no DRC)
                 ssinc = pia * 12.0d0 !translate monthly PIA to yearly ss income 
             endif ! (git < pERA)
          endif    !(giRECSS == 2)
          !====================================================================
          ! earnings test (including Delayed retirement credit)
          !=====================================================================
          if ((gcI_CF .NE. 1) .AND. (git>=pERA) .AND. (git<70)) then !there is earnings test
              if ( ((git<pNRA) .and. (winc>=pETESTbp(1))) .or. ((git>=pNRA) .and. (winc>=pETESTbp(2))) ) then  !before 2000
                  if (git<pNRA) then
                      rtemp = winc - pETESTbp(1) 
                      rtemp = min(ssinc, rtemp / 2.0d0)  !withheld benefit
                  else                 
                      rtemp = winc - pETESTbp(2) 
                      rtemp = min(ssinc, rtemp / 3.0d0)  !withheld benefit
                  endif
                  ssinc = max(0.0d0, ssinc - rtemp)
              endif !if ( ((git<pNRA)
          endif !if (gcI_CF .NE. 1)
      !else  ! no ss in next period (giRECSS_next == 1)
      endif ! if (giHEALTH==4 .AND.

      !3.No Social Security system (no SS tax, no SS benefit).
      !6.Remove SS benefit, but keep SS tax.
      !5.Reduce SS benefit by 20%. SS tax unchanged.
      !10.No Social Security system (no SS tax, no SS benefit)---keep SSDI unchanged
      !11.Reduce SS benefit by 20%. SS tax unchanged.---keep SSDI unchanged
      !12.Remove SS benefit, but keep SS tax;---keep SSDI unchanged
      if ( (gcI_CF==3) .OR. (gcI_CF==6) ) then
          ssinc = 0.0d0
          ssdiinc = 0.0d0
      elseif (gcI_CF == 5) then
          ssinc = 0.8d0 * ssinc
          ssdiinc = 0.8d0 * ssdiinc
      elseif ( (gcI_CF==10) .OR. (gcI_CF==12) ) then
          ssinc = 0.0d0
      elseif (gcI_CF == 11) then
          ssinc = 0.8d0 * ssinc   
      endif     
     
      !positive capital interest income is taxable for income tax, but not for social security tax
      !only earned income, winc, is subject to the social security tax
      
      rtinc = winc + rothinc + max(0.0d0, pr*gA)  !taxable income
      !taxable social security benefit
      if (rtinc <= 0.0d0) then
          MSM_Indiv%a(18,git_pos,i) = 0.0d0 !ss tax
          MSM_Indiv%a(19,git_pos,i) = 0.0d0 !medicare tax
          MSM_Indiv%a(20,git_pos,i) = 0.0d0 !income tax
          MSM_Indiv%a(21,git_pos,i) = ssinc !ss benefit income          
      else
          rtemp = rtinc + (ssinc / 2.0d0)
          if (rtemp <= gr_sstax_1) then
              rtemp2 = 0.0d0 !taxable ss benefit 
          else    
              rtemp2 = min(0.85d0*ssinc, 0.5d0*min(ssinc, rtemp-gr_sstax_1, gr_sstax_2)+0.85d0*max(0.0d0, rtemp-gr_sstax_3))  !taxable ss benefit
          endif   
          rtemp = ssinc - rtemp2  !nontaxable ss benefit part
          rtinc = rtinc + rtemp2  !taxable income, including ss benefit
          !The payroll taxes include the Social Security portion, 6.2% capped at $87, 900, 
          ! and the Medicare tax, 1.45% uncapped.
          !MSM_Indiv%a(18,git_pos,i) = rtax_ss * min(rtinc, rtax_ss_cap) !ss tax
          !MSM_Indiv%a(19,git_pos,i) = rtax_medicare * rtinc  !medicare tax          
          MSM_Indiv%a(18,git_pos,i) = rtax_ss * min(winc, rtax_ss_cap) !ss tax
          MSM_Indiv%a(19,git_pos,i) = rtax_medicare * winc             !medicare tax
          rtemp3 = rtinc - ftaxcode(rtinc) - MSM_Indiv%a(18,git_pos,i) - MSM_Indiv%a(19,git_pos,i)  !federal income tax, excluding SS tax and medicare tax
          rtemp3 = max(0.0d0, rtemp3)
          
          MSM_Indiv%a(20,git_pos,i) = rtemp3 * (rtinc-rtemp2)/rtinc  !federal income tax on wage income, excluding SS tax and medicare tax
          MSM_Indiv%a(21,git_pos,i) = ssinc - rtemp3 * (rtemp2)/rtinc !post-tax ss benefit income, net of taxes on taxable ss income
      endif
      MSM_Indiv%a(22,git_pos,i) = ssdiinc !SSDI & SSI benefit income  !Not taxable
  end subroutine sub_calbudget
  
  !output print. called in master
  subroutine output_simprofile()
     !Print out results
     implicit none
     integer(4) :: i
     real(8) :: rtemp1, rtemp2, rtemp3, rtemp4
     character(100) :: filename, filename_indiv

     !print out moments
     ! and print out individual simulations
     !1-A,2-H,3-AIME,4-c,5-I,6-L,7-Labor,8-Value,9-leisure error
     ! Calculate moments in the output file: 
     !  1 lfpr_mean, 2 lnw_mean, 3 lnw_sd, 4 lnw_fe_mean, 
     !  5 E to U, 6 U to E  7 consumption  8 SS 
     !  9 diff_lfpr_E2G, 10 diff_lfpr_G2B, 11 lnw depreciation after unemployment
     if (gcI_elas_t > 0 .OR. gcI_taxhike_t > 0) then
        ! This is when calculating the elasticity at Age gcI_elas_t
        !do git = N_data_0, N_data_10
        do git = N_0, N_T
           git_pos = git - N_0 + 1 
           !if (MSM_Moments_Sim%d2 == N_moments_type) then
           if (IHasHealthAge > 0) then
              rtemp1 = MSM_Moments_Sim%a(git,7)  !diff_lfpr_E2G
              rtemp2 = MSM_Moments_Sim%a(git,8)  !diff_lfpr_G2B
              rtemp3 = MSM_Moments_Sim%a(git,9)  !diff_lfpr_B2D
           else
              rtemp1 = -999.0d0
              rtemp2 = -999.0d0
              rtemp3 = -999.0d0
           endif
           if (Ipt == 1) then !PT
               if (IHasHealthAge > 0) then
                   rtemp4 = MSM_Moments_Sim%a(git,10)  !pt
               else
                   rtemp4 = MSM_Moments_Sim%a(git,7)  !pt
               endif               
           else
               rtemp4 = -999.0d0
           endif
           
           if (gcI_elas_t > 0) then
               i = gcI_elas_t
           else
               i = gcI_taxhike_t
           endif
           
           write(1100,'(I4,A,I4,A,F12.4,A,F12.4,A,F12.4,A,F12.4,A,F12.4, &
                   & A,F12.4,A,F12.4,A,F12.4,A,F12.4,A,F12.4,   &
                   & A,F12.4,A,F12.4,A,F12.4,A,F12.4,A,F12.4,A,F12.4,A,F12.4,      &
                   & A,F12.4,A,F12.4,A,E40.20,A,E40.20,A,E40.20,A,E40.20,A,E40.20)') i, ',', git, &
                   & ',',MSM_Moments_Sim%a(git,1),   &
                   & ',',MSM_Moments_Sim%a(git,2),   &
                   & ',',MSM_Moments_Sim%a(git,3),   &
                   & ',',MSM_Moments_Sim%a(git,4),   &
                   & ',',MSM_Moments_Sim%a(1,1),     &   !E to U
                   & ',',MSM_Moments_Sim%a(2,1),     &   !U to E
                   & ',',MSM_Moments_Sim%a(git,5),   &   !#5 consumption
                   & ',',MSM_Moments_Sim%a(git,6),   &   !#6 SS
                   & ',',rtemp1, ',',rtemp2,',',rtemp3,',',rtemp4,   &   !diff_lfpr_E2G, diff_lfpr_G2B, diff_lfpr_B2D, lfpr_pt
                   & ',',MSM_Moments_Sim%a(3,1),     &   !wage depreciation after unemployment spell
                   & ',',Mean(MSM_Indiv%a(1,git_pos,:)),   &  !A
                   & ',',Mean(MSM_Indiv%a(2,git_pos,:)),   &  !H
                   & ',',Mean(MSM_Indiv%a(4,git_pos,:)),   &  !c
                   & ',',Mean(MSM_Indiv%a(5,git_pos,:)),   &  !I
                   & ',',Mean(MSM_Indiv%a(7,git_pos,:)),   &  !Labor
                   & ',',Mean(pw(git)*(MSM_Indiv%a(2,git_pos,:)*MSM_Indiv%a(7,git_pos,:))),  &  !Labor income
                   & ',',Mean(MSM_Indiv%a(12,git_pos,:)),   &  !Value of leisure = 0
                   & ',',Mean(MSM_Indiv%a(13,git_pos,:)),   &  !Value of leisure = 1
                   & ',',Mean(MSM_Indiv%a(14,git_pos,:)),   &  !prob of leisure = 0
                   & ',',Mean(MSM_Indiv%a(23,git_pos,:)),   &  !Value of leisure = -1
                   & ',',Mean(MSM_Indiv%a(24,git_pos,:))       !prob of leisure = p
        enddo
        return
     else
        !20 is for moments and 26 is for individual
        if (gcI_CF==88) then
           write(filename, "(A26)") "output/txt_MSM_Moments.txt"
           write(filename_indiv, "(A24)") "output/txt_MSM_Indiv.txt"
        elseif (gcL_CF_IW_decompose .AND. (gcI_CF>=1 .AND. gcI_CF<=12) .AND. (gcI_CF_IW>=1 .AND. gcI_CF_IW<=7)) then
            if (gcI_CF<10) then
                write(filename, "(A25,I1,A2,I1,A4)") "output/txt_MSM_Moments_CF",gcI_CF,"IW",gcI_CF_IW,".txt"
                write(filename_indiv, "(A23,I1,A2,I1,A4)") "output/txt_MSM_Indiv_CF",gcI_CF,"IW",gcI_CF_IW,".txt"
            elseif (gcI_CF<100) then
                write(filename, "(A25,I2,A2,I1,A4)") "output/txt_MSM_Moments_CF",gcI_CF,"IW",gcI_CF_IW,".txt"
                write(filename_indiv, "(A23,I2,A2,I1,A4)") "output/txt_MSM_Indiv_CF",gcI_CF,"IW",gcI_CF_IW,".txt"
            endif
        elseif (gcI_CF<10) then
           write(filename, "(A25,I1,A4)") "output/txt_MSM_Moments_CF",gcI_CF,".txt"
           write(filename_indiv, "(A23,I1,A4)") "output/txt_MSM_Indiv_CF",gcI_CF,".txt"
        elseif (gcI_CF<100) then
           write(filename, "(A25,I2,A4)") "output/txt_MSM_Moments_CF",gcI_CF,".txt"
           write(filename_indiv, "(A23,I2,A4)") "output/txt_MSM_Indiv_CF",gcI_CF,".txt"
        elseif (gcI_CF<1000) then
           write(filename, "(A25,I3,A4)") "output/txt_MSM_Moments_CF",gcI_CF,".txt"
           write(filename_indiv, "(A23,I3,A4)") "output/txt_MSM_Indiv_CF",gcI_CF,".txt"   
        endif
        open(20,file=trim(cwd)//trim(filename))
        open(26,file=trim(cwd_indiv)//trim(filename_indiv))

        !===========================================================================================================================
        ! Calculate moments: 1 lfpr_mean, 2 lnw_mean, 3 lnw_sd, 4 lnw_fe_mean, 5 E to U, 6 U to E
        !  7 consumption  8 diff_lfpr
        !===========================================================================================================================
        !do git = N_data_0, N_data_10
        do git = N_0, N_T
           if (IHasHealthAge > 0) then
              rtemp1 = MSM_Moments_Sim%a(git,7)  !diff_lfpr_E2G
              rtemp2 = MSM_Moments_Sim%a(git,8)  !diff_lfpr_G2B
              rtemp3 = MSM_Moments_Sim%a(git,9)  !diff_lfpr_B2D
           else
              rtemp1 = -999.0d0
              rtemp2 = -999.0d0
              rtemp3 = -999.0d0
           endif
           if (Ipt == 1) then !PT
               if (IHasHealthAge > 0) then
                   rtemp4 = MSM_Moments_Sim%a(git,10)  !pt
               else
                   rtemp4 = MSM_Moments_Sim%a(git,7)  !pt
               endif               
           else
               rtemp4 = -999.0d0
           endif
            
           write(20,'(I4,A,F12.4,A,F12.4,A,F12.4,A,F12.4,A,F12.4,     & 
                 &       A,F12.4,A,F12.4,A,F12.4,A,F12.4,A,F12.4,A,F12.4,A,F12.4,A,F12.4)') git,  &
                   & ',',MSM_Moments_Sim%a(git,1),   &
                   & ',',MSM_Moments_Sim%a(git,2),   &
                   & ',',MSM_Moments_Sim%a(git,3),   &
                   & ',',MSM_Moments_Sim%a(git,4),   &
                   & ',',MSM_Moments_Sim%a(1,1),     &  !E to U
                   & ',',MSM_Moments_Sim%a(2,1),     &  !U to E
                   & ',',MSM_Moments_Sim%a(git,5),   &  !consumption
                   & ',',MSM_Moments_Sim%a(git,6),   &  !#6 SS
                   & ',',rtemp1,',',rtemp2,',',rtemp3,',',rtemp4,    &  !diff_lfpr_E2G, diff_lfpr_G2B, diff_lfpr_B2D, pt
                   & ',',MSM_Moments_Sim%a(3,1)         !wage depreciation after unemployment spell
        enddo
        close(20)
     endif   !if (gcI_elas_t > 0) then     
     
     do i = 1, N_Sim
        do git = N_0, N_T
           git_pos = git - N_0 + 1 
           !1-A,2-H,3-AIME,4-c,5-I,6-L,7-Labor,8-Value,9-leisure error,10-XI error, 11-HEALTH shock
           write(26,'(I10,A,I10,A, &
                    & E100.50,A,E100.50,A,E100.50,A,E100.50,A,E100.50,  &
                    & A,E100.50,A,E100.50,A,E100.50,A,E100.50,A,E100.50,    &
                    & A,E100.50,A,E100.50,A,E100.50,A,E100.50,A,E100.50,A,E100.50,  &
                    & A,I10,A,I10,A,I10,A,I10,A,I10,A,I10,A,I10,  &
                    & A,E100.50,A,E100.50,A,E100.50,A,E100.50,A,E100.50,A,E100.50,A,E100.50,A,E20.10)')  &
                                 & i,',',git, &
                                 & ',',MSM_Indiv%a(1,git_pos,i), &
                                 & ',',MSM_Indiv%a(2,git_pos,i), &
                                 & ',',MSM_Indiv%a(3,git_pos,i), &
                                 & ',',MSM_Indiv%a(4,git_pos,i), &
                                 & ',',MSM_Indiv%a(5,git_pos,i), &
                                 & ',',MSM_Indiv%a(6,git_pos,i), &    !Leisure, 0/1 or 0/0.5/1
                                 & ',',MSM_Indiv%a(7,git_pos,i), &
                                 & ',',MSM_Indiv%a(8,git_pos,i), &    !value, including leisure
                                 & ',',MSM_Indiv%a(9,git_pos,i), &    !leisure error --- standard normal
                                 & ',',MSM_Indiv%a(10,git_pos,i), &
                                 & ',',MSM_Indiv%a(11,git_pos,i), &
                                 & ',',MSM_Indiv%a(12,git_pos,i), &   !V0
                                 & ',',MSM_Indiv%a(13,git_pos,i), &   !V1
                                 & ',',MSM_Indiv%a(14,git_pos,i), &   !prob of leisure=0
                                 & ',',MSM_Indiv%a(15,git_pos,i), &   !FS transition error
                                 & ',',MSM_Indiv%a(16,git_pos,i), &   !FS spouse actual inc
                                 & ',',MSM_Indiv_Irecss%a(git_pos,i), &  !ss status
                                 & ',',MSM_Indiv_Irecss%a(min(N_T-N_0+1,git_pos+1),i), &  !ss policy
                                 & ',',MSM_Indiv_Ihealth%a(git_pos,i), &    !health status 
                                 & ',',MSM_Indiv_IFS%a(git_pos,i),  &        !Family status
                                 & ',',MSM_Indiv_Ihet%a(i,1),     &     !het type
                                 & ',',MSM_Indiv_Ihet%a(i,2),     &     !observed first year
                                 & ',',MSM_Indiv_Ihet%a(i,3),     &     !observed last year
                                 & ',',MSM_Indiv%a(17,git_pos,i),   &     !NEW: 17-produced output
                                 & ',',MSM_Indiv%a(18,git_pos,i),   &     !18-ss tax
                                 & ',',MSM_Indiv%a(19,git_pos,i),   &     !19-medicare tax
                                 & ',',MSM_Indiv%a(20,git_pos,i),   &     !20-income tax, excluding ss tax and medicare tax
                                 & ',',MSM_Indiv%a(21,git_pos,i),   &     !21-received ss benefit
                                 & ',',MSM_Indiv%a(22,git_pos,i),   &     !22-ssdiinc
                                 & ',',MSM_Indiv%a(23,git_pos,i),   &     !23-Vp
                                 & ',',MSM_Indiv%a(24,git_pos,i)          !24-prob of leisure=p
                                 

        enddo !git
     enddo  !i
     close(26)
     
  end subroutine output_simprofile

  
  !========================================================================================================
  ! calculate the alternative situation
  ! Simulate lifecycle profiles all by interpolation, called by all nodes (master, head slaves and slaves)
  !========================================================================================================
  subroutine SMM_SimulateProfiles_Alt()  !!!!!!NEED TO FIX FOR Ipt=1
    implicit none
    integer(4) :: i, j, ix(3), ipos, ipos2
    real(8) :: rtemp, rtemp2, rnext(2), rcx3(3), rcx2(2), ev, rAfix, rAIMEfix
    
    MSM_Indiv_Alt%d1 = MSM_Indiv%d1
    MSM_Indiv_Alt%d2 = N_T_grid
    MSM_Indiv_Alt%d3 = N_Sim
    ALLOCATE(MSM_Indiv_Alt%a(MSM_Indiv_Alt%d1,MSM_Indiv_Alt%d2,MSM_Indiv_Alt%d3))
    MSM_Indiv_Alt%a = 0.0d0
    gI_Indiv_Alt_size_mpi = MSM_Indiv_Alt%d1*MSM_Indiv_Alt%d2

    MSM_Indiv_Irecss_Alt%d1 = 2    !leisure 1 and 0
    MSM_Indiv_Irecss_Alt%d2 = N_T_grid
    MSM_Indiv_Irecss_Alt%d3 = N_Sim
    ALLOCATE(MSM_Indiv_Irecss_Alt%a(MSM_Indiv_Irecss_Alt%d1, MSM_Indiv_Irecss_Alt%d2, MSM_Indiv_Irecss_Alt%d3))
    MSM_Indiv_Irecss_Alt%a = 0
        
    !MSM_Indiv_Alt
    !1-c,2-I,3-L,4-A(end of git),5-AIME prime,6-EV !Leisure=1
    !7-c,8-I,9-L,10-A(end of git),11-AIME prime,12-EV,  !Leisure=0
    !13-EV(A&H no change),14-EV(A&AIME no change),15-EV(A no change)  !Leisure=0
    !16-part 1: wage income
    !17-part 2: from increasing AIME only (13-6)
    !18-part 3: from increasing H only (14-6)
    !19-from increasing AIME and H: 15-6
    !20-v1(Leisure=1), 21-v1(Leisure=1, A+wage income), 22-value of increasing A

        
    !do i = 1, N_Sim
    do i = gI_Indiv_Tasks_mpi%a(myid_mpi+1,1), gI_Indiv_Tasks_mpi%a(myid_mpi+1,2) 
              
       do git = N_0, N_T
          git_pos = git - N_0 + 1 
          call xPosadjust(git, MSM_Indiv%a(1:3,git_pos,i), ix)
          gA = MSM_Indiv%a(1,git_pos,i)
          gH = MSM_Indiv%a(2,git_pos,i)
          gAIME = MSM_Indiv%a(3,git_pos,i)
          giA = ix(1)
          giH = ix(2)
          giAIME = ix(3)
          giRECSS = MSM_Indiv_Irecss%a(git_pos,i) + 1   !index, 1-not RECSS, 2-RECSS
          giHEALTH = MSM_Indiv_Ihealth%a(git_pos,i)      !health status, 1-excellent, 2-good, 3-bad 4-disabled
          giFS = MSM_Indiv_IFS%a(git_pos,i)   !Family status: 1-3
        
          !=========================================
          !the optimal solution in the simulation
          !=========================================
          ipos = 1 - NINT(MSM_Indiv%a(6,git_pos,i))  !0 if leisure==1, 1 if leisure==0
          ipos2 = ipos * 6
          if (git < N_T) then
              MSM_Indiv_Irecss_Alt%a(1+ipos, git_pos+1,i) = MSM_Indiv_Irecss%a(git_pos+1,i)  !0 or 1
              giRECSS_next = 1 + MSM_Indiv_Irecss_Alt%a(1+ipos, git_pos+1,i)
          else
              giRECSS_next = 2 !always receive SS at age N_T
          endif
          MSM_Indiv_Alt%a(1+ipos2:3+ipos2,git_pos,i) = MSM_Indiv%a(4:6,git_pos,i)  !c,I,L
          !The next period state variables
          rtemp2 = MSM_Indiv_Alt%a(1+ipos2,git_pos,i)  !c
          rtemp = pw(git)*gH*(1.0d0-MSM_Indiv_Alt%a(2+ipos2,git_pos,i))*(1.0d0-MSM_Indiv_Alt%a(3+ipos2,git_pos,i))  !Labor income
          call A_AIME_prime(rtemp2, rtemp, rnext)  !given git, gA, gAIME, giRECSS, giRECSS_next (current policy) 
          MSM_Indiv_Alt%a(6+ipos2,git_pos,i) = cmaxee_ev(MSM_Indiv_Alt%a(2+ipos2:3+ipos2,git_pos,i), rnext)  !EV
          MSM_Indiv_Alt%a(4+ipos2,git_pos,i) = rnext(1)  !A prime
          MSM_Indiv_Alt%a(5+ipos2,git_pos,i) = rnext(2)  !AIME prime
                    
          !=====================================================
          !now we calculate the other option, leisure = j
          ! and assume the SSA is the same as the previous case
          !=====================================================
          j = 1 - NINT(MSM_Indiv%a(6,git_pos,i)) !leisure is j, the opposite of simulated choice
          ipos = 1 - j   !0 if j==1, 1 if j==0
          if (git < N_T) MSM_Indiv_Irecss_Alt%a(1+ipos, git_pos+1,i) = MSM_Indiv_Irecss%a(git_pos+1,i)  !0 or 1
          ipos2 = ipos * 6
          MSM_Indiv_Alt%a(3+ipos2,git_pos,i) = dble(j)  !leisure
          !deciding optimal choice now
          call sub_ForwardInt(MSM_Indiv_Alt%a(1+ipos2:3+ipos2,git_pos,i),rnext,ev)
          MSM_Indiv_Alt%a(6+ipos2,git_pos,i) = cmaxee_ev(MSM_Indiv_Alt%a(2+ipos2:3+ipos2,git_pos,i), rnext)  !EV
          MSM_Indiv_Alt%a(4+ipos2,git_pos,i) = rnext(1)  !A prime
          MSM_Indiv_Alt%a(5+ipos2,git_pos,i) = rnext(2)  !AIME prime
          
          !=====================================================
          !now we calculate the hyperthetical case
          !=====================================================
          if (git < N_T) then
              giRECSS_next = MSM_Indiv_Irecss_Alt%a(1,git_pos+1,i) + 1  !use Leisure=1 as the base case
          else
              giRECSS_next = 2
          endif
          rAfix = MSM_Indiv_Alt%a(4,git_pos,i)      !A prime for Leisure=1
          rAIMEfix = MSM_Indiv_Alt%a(5,git_pos,i)   !AIME prime for Leisure=1
          
          !in increasing AIME only, keeping A&H unchanged
          rcx2(1) = 0.0d0   !I
          rcx2(2) = 1.0d0  !assuming Leisure=1 so no investment in H 
          rnext(1) = rAfix
          rnext(2) = MSM_Indiv_Alt%a(11,git_pos,i)  !AIME prime for Leisure=0
          MSM_Indiv_Alt%a(13,git_pos,i) = cmaxee_ev(rcx2, rnext)  !13-EV(A&H no change)
              
          !in increasing H only, keeping A&AIME unchanged
          rcx2(1) = MSM_Indiv_Alt%a(8,git_pos,i)   !I for Leisure=0
          rcx2(2) = MSM_Indiv_Alt%a(9,git_pos,i)   !L for Leisure==0
          rnext(1) = rAfix
          rnext(2) = rAIMEfix  !AIME no change, as in Leisure=1
          MSM_Indiv_Alt%a(14,git_pos,i) = cmaxee_ev(rcx2, rnext)  !14-EV(A&AIME no change)
          
          !in increasing AIME and H, keeping A unchanged
          rcx2(1) = MSM_Indiv_Alt%a(8,git_pos,i)   !I for Leisure=0
          rcx2(2) = MSM_Indiv_Alt%a(9,git_pos,i)   !L for Leisure==0
          rnext(1) = rAfix
          rnext(2) = MSM_Indiv_Alt%a(11,git_pos,i)  !AIME for Leisure=0
          MSM_Indiv_Alt%a(15,git_pos,i) = cmaxee_ev(rcx2, rnext)  !15-EV(A no change)

          !20-v1(Leisure=1), 21-v0(Leisure=0, only from A), 22-value of increasing A
          MSM_Indiv_Alt%a(20,git_pos,i) = linear_inte_tetrahedra3(sgridA%a(git_pos,ix(1):(ix(1)+1)), &
                         &  sgridAIME%a(git_pos,ix(3):(ix(3)+1)), sgridH%a(git_pos,ix(2):(ix(2)+1)), &
                         &  dvV(1)%a(ix(1):(ix(1)+1),ix(3):(ix(3)+1),giHEALTH,giRECSS,giFS,ix(2):(ix(2)+1),git_pos),  &
                         &  (/gA,gAIME,gH/) ) !Value except leisure part for Leisure=1
          !labor income for Leisure=0
          rtemp = pw(git)*gH*(1.0d0-MSM_Indiv_Alt%a(8,git_pos,i))  !Labor income
          rcx3(1) = gA + ftaxcode(rtemp)    !A
          rcx3(2) = gH                      !H
          rcx3(3) = gAIME                   !AIME
          call xPosadjust(git, rcx3, ix)
          MSM_Indiv_Alt%a(21,git_pos,i) = linear_inte_tetrahedra3(sgridA%a(git_pos,ix(1):(ix(1)+1)), &
                         &  sgridAIME%a(git_pos,ix(3):(ix(3)+1)), sgridH%a(git_pos,ix(2):(ix(2)+1)), &
                         &  dvV(1)%a(ix(1):(ix(1)+1),ix(3):(ix(3)+1),giHEALTH,giRECSS,giFS,ix(2):(ix(2)+1),git_pos),  &
                         &  (/rcx3(1), rcx3(3), rcx3(2)/) ) !Value except leisure part for Leisure==1 with higher gA
          
          
          !Now calculate u'_c, the Leisure=1 as the base
          gI_p_L = NINT(MSM_Indiv_Alt%a(3,git_pos,i))          
          rtemp = gr_psi(git, giFS, gI_p_L) * (MSM_Indiv_Alt%a(1,git_pos,i)**(-petac))   !leisure=1
          
          rtemp = pbeta / rtemp  !\tilde{beta}
          !11-part 1: wage income
          !12-part 2: from increasing AIME only
          !13-part 3: from increasing H only
          !14-total change in EV: 8-4
          MSM_Indiv_Alt%a(16,git_pos,i) = pw(git)*gH*(1.0d0-MSM_Indiv_Alt%a(8,git_pos,i))  !1-income effect (wage income)
          MSM_Indiv_Alt%a(17,git_pos,i) = MSM_Indiv_Alt%a(13,git_pos,i) - MSM_Indiv_Alt%a(6,git_pos,i)
          MSM_Indiv_Alt%a(18,git_pos,i) = MSM_Indiv_Alt%a(14,git_pos,i) - MSM_Indiv_Alt%a(6,git_pos,i)
          MSM_Indiv_Alt%a(19,git_pos,i) = MSM_Indiv_Alt%a(15,git_pos,i) - MSM_Indiv_Alt%a(6,git_pos,i)
          MSM_Indiv_Alt%a(17:19,git_pos,i) = rtemp * MSM_Indiv_Alt%a(17:19,git_pos,i)
          MSM_Indiv_Alt%a(22,git_pos,i) = MSM_Indiv_Alt%a(21,git_pos,i) - MSM_Indiv_Alt%a(20,git_pos,i)
          MSM_Indiv_Alt%a(22,git_pos,i) = rtemp * MSM_Indiv_Alt%a(22,git_pos,i)
          
      enddo !do git = N_0, N_T
    enddo !i
    
    !to print out simulated individual results
    i = min(gI_Indiv_Tasks_mpi%a(myid_mpi+1,1), N_Sim) 
    call MPI_Gatherv(MSM_Indiv_Alt%a(1,1,i), gI_Indiv_Alt_size_mpi*gI_Indiv_Tasks_mpi%a(myid_mpi+1,3), MPI_DOUBLE_PRECISION,   & 
                  &  MSM_Indiv_Alt%a(1,1,1), gI_Indiv_Alt_size_mpi*gI_Indiv_Tasks_mpi%a(:,3), &
                  &  gI_Indiv_Alt_size_mpi*gI_Indiv_Tasks_mpi%a(:,4),      &
                  &  MPI_DOUBLE_PRECISION, master_mpi, MPI_COMM_WORLD, ipos)
    call MPI_Gatherv(MSM_Indiv_Irecss_Alt%a(1,1,i), 2*N_T_grid*gI_Indiv_Tasks_mpi%a(myid_mpi+1,3), MPI_INTEGER,   & 
                  &  MSM_Indiv_Irecss_Alt%a(1,1,1), 2*N_T_grid*gI_Indiv_Tasks_mpi%a(:,3),  &
                  &  2*N_T_grid*gI_Indiv_Tasks_mpi%a(:,4),     &
                  &  MPI_INTEGER, master_mpi, MPI_COMM_WORLD, ipos)

    if (myid_mpi == master_mpi) call sub_Indiv_Alt_Output()
    
    DEALLOCATE(MSM_Indiv_Alt%a)
    DEALLOCATE(MSM_Indiv_Irecss_Alt%a)
 end subroutine SMM_SimulateProfiles_Alt
  

   !called in master only. run in master only.
 subroutine sub_Indiv_Alt_Output()
    implicit none
    integer(4) :: i

    open(260,file=trim(cwd_indiv)//"output/txt_MSM_Indiv_Alt.txt")
        
    do i = 1, N_Sim
        do git = N_0, N_T
           git_pos = git - N_0 + 1 
           write(260,'(I10,A,I10,A, &
                    & E100.50,A,E100.50,A,E100.50,A,E100.50,A,E100.50,  &
                    & A,E100.50,A,E100.50,A,E100.50,A,E100.50,A,E100.50,    &
                    & A,E100.50,A,E100.50,A,E100.50,A,E100.50,A,E100.50,    &
                    & A,E100.50,A,E100.50,A,E100.50,                        &
                    & A,E100.50,A,E100.50,A,E100.50,A,E100.50,A,E100.50,A,E100.50)')  &
                                 & i,',',git, &
                                 & ',',MSM_Indiv_Alt%a(1,git_pos,i), &
                                 & ',',MSM_Indiv_Alt%a(2,git_pos,i), &
                                 & ',',MSM_Indiv_Alt%a(3,git_pos,i), &
                                 & ',',MSM_Indiv_Alt%a(4,git_pos,i), &
                                 & ',',MSM_Indiv_Alt%a(5,git_pos,i), &
                                 & ',',MSM_Indiv_Alt%a(6,git_pos,i), &
                                 & ',',MSM_Indiv_Alt%a(7,git_pos,i), &
                                 & ',',MSM_Indiv_Alt%a(8,git_pos,i), &
                                 & ',',MSM_Indiv_Alt%a(9,git_pos,i), &
                                 & ',',MSM_Indiv_Alt%a(10,git_pos,i), &
                                 & ',',MSM_Indiv_Alt%a(11,git_pos,i), &
                                 & ',',MSM_Indiv_Alt%a(12,git_pos,i), &
                                 & ',',MSM_Indiv_Alt%a(13,git_pos,i), &
                                 & ',',MSM_Indiv_Alt%a(14,git_pos,i), &
                                 & ',',MSM_Indiv_Alt%a(15,git_pos,i), &
                                 & ',',MSM_Indiv_Alt%a(16,git_pos,i), &
                                 & ',',MSM_Indiv_Alt%a(17,git_pos,i), &
                                 & ',',MSM_Indiv_Alt%a(18,git_pos,i), &
                                 & ',',MSM_Indiv_Alt%a(19,git_pos,i), &
                                 & ',',MSM_Indiv_Alt%a(20,git_pos,i), &
                                 & ',',MSM_Indiv_Alt%a(21,git_pos,i), &
                                 & ',',MSM_Indiv_Alt%a(22,git_pos,i), &
                                 & ',',MSM_Indiv_Alt%a(23,git_pos,i),   &     !23-Vp
                                 & ',',MSM_Indiv_Alt%a(24,git_pos,i)          !24-prob of leisure=p
        enddo !git
     enddo  !i
     close(260)
  end subroutine sub_Indiv_Alt_Output

   
END MODULE DISTANCE
